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CHAPTER  1: 

Introduction 


As  the  Department  of  Defense  (DOD)  continues  its  development  and  implementation  of 
a  new  category  of  weaponry,  directed  energy  weapons  (DEW),  other  nations  are  expected 
to  follow  suit  [1],  The  design  method  known  as  transformation  optics  (TO),  developed  by 
Pendry,  Schurig,  and  Smith  [2],  provides  a  highly  useful  design  tool  for  determining  the 
electromagnetic  material  properties  required  in  order  to  redirect  incoming  electromagnetic 
radiation  along  a  more  desirable  path.  This  design  method  has  proven  effective  in  the 
fabrication  of  real  world  structures  using  metamaterials  —  materials  engineered  to 
possess  properties  which  are  not  found  in  nature  [3].  Although  both  powerful  and 
elegant  in  its  formulation,  TO  relies  upon  a  purely  linear  formulation  of  the  response 
of  materials  to  electric  and  magnetic  fields.  It  has  been  known  since  at  least  1941  [4] 
that  sufficiently  large  electric  and  magnetic  field  amplitudes  induce  nonlinear  (i.e.,  second 
and  higher  order)  material  responses.  As  these  high  field  amplitudes  are  naturally  expected 
in  an  environment  in  which  the  concern  is  the  countering  of  certain  classes  of  DEW, 
such  as  high  power  lasers  and  microwave  devices,  a  question  of  the  effectiveness  of 
TO  derived  redirection  structures  for  counter-directed  energy  weapon  (CDEW) 
applications  arises.  We  attempt  here  to  address  this  question  of  the  utility  of  the 
transformation  optics  approach  in  a  CDEW  environment. 


1.1  Directed  Energy  Weapons 

A  directed  energy  weapon  is  a  weapon  which  delivers  electromagnetic  energy,  rather  than 
a  traditional  projectile,  to  its  target.  Typical  examples  of  directed  energy  weapons  include 
high  power  microwave  (HPM)  or  high  power  radio  frequency  (HPRF)  devices,  designed  to 
induce  large  voltages  and  currents  in  adversary  electronic  systems  in  order  to  render  those 
systems  inoperative,  and  high  energy  lasers  (HEL),  which  ablate  target  material  through 
the  delivery  of  a  focused  high  intensity  beam  of  electromagnetic  radiation  onto  a  small 
area  of  the  target.  These  weapon  systems  have  been  active  areas  of  Navy  research 
since  the  1960s,  and  are  the  subject  of  ongoing  development  [1],  [5],  As  research  into 
these  devices  proceeds,  concern  has  grown  regarding  the  ability  of  potential  adversaries 
to  field  similar 
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weaponry,  particularly  within  the  realm  of  anti  satellite  applications  [6]. 


1.2  Counter-Directed  Energy  Weapons 

In  response  to  this  concern,  the  Office  of  Naval  Research  (ONR),  in  conjunction  with  the 
Naval  Postgraduate  School  (NPS),  United  States  Naval  Academy  (USNA),  and  Naval  Re¬ 
search  Laboratory  (NRL)  is  "investigating  basic  research  topics  in  counter  threats  from 
directed  energy  weapons  systems,  such  as  high-power  lasers  or  microwaves"  [7].  These 
research  topics  include,  from  [1], 

a.  Advanced  materials  including  nano-  and/or  nonlinear  materials  for  enhanced 
HEL  protection  of  sensors,  optics,  airframe,  etc. 

b.  Metamaterial  structures  for  the  control  and  mitigation  of  HEL  and  HPRF 
irradiation. 

c.  Techniques  for  HEL  mitigation  such  as  use  of  plasmas  and  obscurants. 

d.  HEL  protection  by  degrading  atmospheric  transmission  (e.g.  thermal 
blooming,  scattering,  absorption  aids,  and  turbulence). 

e.  Modeling  and  sensing  of  off-axis  detection  and  source  geo-location. 

f.  Novel  instrumentation  for  detection  of  HEL  and  HPRF  irradiation. 

g.  Active/passive  circuit  protection  and  limiters  for  HPRF 

h.  Modeling  of  HPRF  and  HEL  effects  to  materials,  electronics  and  sensors  as 
applied  to  CDEW  objectives. 


1.2.1  Transformation  Optics 

The  field  of  transformation  optics,  founded  in  [2]  by  Pendry,  Schurig,  and  Smith,  uses  the 
form  invariance  of  Maxwell’s  equations,  the  fundamental  equations  of  electromagnetism, 
as  well  as  the  transformation  properties  of  various  electromagnetic  fields  and  material  prop¬ 
erties  under  coordinate  transformations,  in  order  to  determine  the  material  properties  nec¬ 
essary  to  force  electromagnetic  fields  to  behave  in  some  desired  fashion.  A  design  tool  for 
directing  electromagnetic  radiation  along  a  specified  path,  transformation  optics’  potential 
for  CDEW  applications  was  quickly  recognized.  Just  as  quickly,  the  material  parameters 
necessary  to  realize  a  working  electromagnetic  redirection  structure  (as  seen  below),  were 
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discovered  to  be  unlike  those  found  in  nature.  Specifically,  negative  electric  and  magnetic 
susceptibility  values  are  commonly  required  [2],  [3],  [8],  [9],  [10].  In  addition,  several 
applications  require  unrealistically  diverging  permittivity  and  permeability  as  an  outer  or 
inner  edge  of  the  structure  is  approached,  as  well  as  continuously  (spatially)  varying  permit¬ 
tivity  and  permeability  [8],  [9].  Figure  1.1  illustrates  both  ideal  and  approximate  material 
parameters  required  to  achieve  the  cylindrical  redirection  structure  of  [9],  as  well  as  the 
simulated  performance  of  the  same  structure  with  ideal  material  parameters. 


Figure  1.1:  Left:  Required  ideal  material  properties  for  a  cylindrical  cloak  of  inner  radius  R i 
and  outer  radius  of  R2  =  2R\,  as  well  as  design  compromises  in  order  to  make  the  structure 
realizeable.  Right:  Full-wave  simulation  of  the  performance  of  the  idealized  material  parameters. 
Note  that  the  waveform  is  undisturbed  upon  exiting  the  structure.  Grey  lines  indicate  the  local 
direction  of  power  flow.  Source:  [11]. 


While  approximations  and  informed  design  compromises  may  effectively  address  the  issue 
of  diverging  or  continuously  varying  material  properties  [3],  [12],  the  problem  of  material 
properties  which  are  not  observed  in  nature  may  be  resolved  through  the  use  of  metamate¬ 
rials  [12],  [13].  Of  particular  note,  however  is  the  completely  linear  constitutive  relation 


Pi  =  e0XijEj  (1.1) 

Mi  =  xiiHj  (1.2) 

employed  in  TO,  where  P  is  the  polarization,  x'e  are  the  components  of  the  electric  sus- 
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ceptibility,  E  is  the  electric  field,  M  is  the  magnetization,  x‘m  are  the  components  of  the 
magnetic  susceptibility,  and  H  is  the  magnetic  field,  and  the  Einstein  summation  conven¬ 
tion  is  in  use.  A  treatment  of  the  theory  of  transformation  optics  is  provided  in  Chapter  3  of 
this  work. 

1.2.2  Metamaterials 

A  metamaterial  is  a  structure  which  is  designed  to  possess  aggregate  material 
properties  which  have  not  yet  been  observed  in  nature.  Within  the  context  of 
electromagnetic  interac-tions,  these  novel  material  responses  are  typically  designed 
through  the  addition  of  struc-tural  features  which  are  smaller  than  the  wavelength  of 
radiation  to  be  scattered.  Specific  to  the  negative  permittivity  and  permeability  often 
required  for  implementation  of  TO  applications,  an  array  of  split-ring  resonators,  structures 
featuring  two  conducting,  concentric,  incomplete  rings,  is  commonly  used.  In  2009,  a 
cylindrical  cloak  operating  at  microwave  frequencies,  designed  through  TO  and  employing 
an  array  of  split-ring  resonators  among  other  design  features,  was  successfully 
demonstrated  [3].  See  Figure  1.2  for  an  example  of  this  design. 


4.7  4.8  4.9  5  5.1  5.2  5.3  5.4  5.5 


Frequency  (GHz) 

Figure  1.2:  Left:  Diagram  of  a  split-ring  resonator  alongside  its  resonance  curve  (source:  [14]). 
Right:  Cylindrical  microwave  frequency  cloak  (source:  [3]),  with  diagrams  of  split-ring  resonator 
design  based  on  radial  position,  and  diagram  of  idealized  material  parameters. 
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Metamaterials  may  therefore  be  viewed  as  a  promising  set  of  structures  with  which  we  may 
realize  a  CDEW  focused  electromagnetic  redirection  device  that  has  been  designed  through 
transformation  optics. 


1.3  Nonlinear  Optics 

As  mentioned  in  the  previous  section,  the  TO  design  approach  assumes  a  linear  set  of 
constitutive  relations  as  given  by  (1.1)  and  (1.2).  High  amplitude  electric  and  magnetic 
fields,  however,  are  known  to  induce  second  order  and  higher  effects  typically  modeled  by, 
as  an  extension  of  (1.1), 


Pl  =  £0  [xijEj  +  xijkEjEk  +  xijk,EjEkE,  + ...]  (1.3) 

where  Xe  are  the  components  of  the  second-order  electric  susceptibility  and  Xe  arc  the 
components  of  the  third-order  electric  susceptibility  [15].  Second-order  effects  include  fre¬ 
quency  doubling,  while  third  order  effects  include  the  Kerr  effect,  wherein  refractive  index 
varies  with  field  strength  [15].  The  linear  constitutive  relations  (1.1)  and  (1.2)  do  not  allow 
for  explicit  realization  of  these  observed  phenomena.  Therefore,  the  TO  approach  does  not 
take  nonlinear  optical  effects  into  account,  and  thus  structures  designed  through  TO  may 
be  expected  to  stray  from  design  expectations  as  field  strengths  rise  high  enough  for  non¬ 
linear  effects  to  become  non-negligible.  As  a  potential  DEW  attack  upon  a  TO  designed 
redirection  structure  may  involve  electric  and  magnetic  fields  of  sufficient  intensity  to  trig¬ 
ger  significant  nonlinear  effects,  the  question  of  applicability  of  the  TO  design  method  for 
CDEW  applications  arises.  This  question  is  the  focus  of  the  following  chapters. 
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CHAPTER  2: 

Classical  Electrodynamics  in  Matter 


As  TO  describes  a  design  method  through  which  we  may  control  the  behavior  of  elec¬ 
tromagnetic  fields  in  matter,  a  description  of  the  physical  phenomenon  of 
electromagnetic  waves  is  in  order.  The  objective  of  the  following  chapter  will  be  the 
description  of  gen-eral  electromagnetic  phenomena  via  the  Maxwell  equations,  followed 
by  the  application  of  these  concepts  to  electromagnetic  waves  in  free  space  as  well  as  in 
media. 


2.1  Maxwell’s  Equations 

Maxwell’s  equations,  the  governing  differential  equations  for  all  classical  electrodynamic 
phenomena,  are  expressed  as: 


W-E  =  — 

£o 

(2.1) 

„  -  dB 

V  xE  =  — T- 
dt 

(2.2) 

V-B  =  0 

(2.3) 

dE 

V  x  B  —  \AqJ  +  Sofia— 

(2.4) 

where  E  is  the  electric  field,  B  the  magnetic  flux  density,  p  the  charge  density,  J  the  cur¬ 
rent  density,  £0  the  permittivity  of  vacuum,  and  fl0  the  permeability  of  vacuum.  Note  that 
overall  charge  conservation  is  consistent  with  the  Maxwell  equations  as  written:  Taking  the 
divergence  of  the  Ampere-Maxwell  law  (2.4), 


V-  (^VxfiJ  =  p0V ■  /  +  e0pov ■ 


Since  the  divergence  of  a  curl  is  always  equal  to  zero,  and  the  spatial  and  time  derivatives 
are  independent  of  each  other, 
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diV -E 

0  =  jUqV  ■  /  +  £0iU0  V  d 


(2.6) 


Inserting  Gauss’  law  (2.1)  into  the  above  relation, 


£o 


dt 


-V-/ 


dp 

dt 


— V-7 


(2.7) 


(2.8) 


2.2  Maxwell’s  Equations  in  Matter 

While  the  Maxwell  equations  (2.1)  -  (2.4)  provide  a  valid  description  of  all  electromagnetic 
phenomena,  they  are  somewhat  cumbersome  when  considering  these  fields  in  matter,  where 
charge  may  be  bound  up  in  a  material’s  constituent  atoms.  A  more  convenient  approach 
is  to  first  draw  a  distinction  between  pj,  the  free  charge  density,  and  p/;,  the  bound  charge 
density.  Defining  the  polarization  as 


p=l%Pl  <2'9> 

where  pj  is  the  electric  dipole  moment  of  a  particular  atom  or  molecule  of  interest.  The 
bound  charge  density  is  defined  as 


Ph  =  ~ V  ■  P 


(2.10) 


If  charge  is  to  be  conserved,  then 


dpb 

dt 


-V-?p 


(2.11) 
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where  Jp  is  the  polarization  current  associated  with  the  motion  of  bound  charges.  The 
above  relationship  thus  requires 


Jn 


dP 

~di 


yielding 


dpb 

dt 


=  -V-7p 


(2.12) 


(2.13) 


The  free  charge  density  pj  is  defined  as  any  charge  density  which  is  not  part  of  the  bound 
charge  density,  giving 


pf  +  pb  =  p  (2.14) 

Similarly,  the  magnetization  M  is  defined 

M=^m  (2.15) 

where  m,  is  the  magnetic  moment  of  an  atom  or  molecule  of  interest.  The  bound  current 
density  Jb  is 

Jb  =  VxM  (2.16) 

The  free  current  density  Jj-  is  the  current  density  unaccounted  for  by  the  bound  or  polariza¬ 
tion  current  densities.  That  is, 


J  —  J f  +  Jp  +  Jb 


(2.17) 
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Inserting  (2.14)  into  (2.1), 


S/E  =  Pf  +  Ph 
£o 


(2.18) 


S/E 


Pf  V-P 

£()  £o 


(2.19) 


V-(e0  E  +  Pj=pf  (2.20) 

The  quantity  £qE  +  P  is  the  electric  displacement  field,  denoted  D.  Arriving  at  Gauss’  Law 
for  macroscopic  media, 


S7  D  =  pf 


Inserting  (2.17)  into  (2.4), 


„  ,  dE 

V  X  B  —  Hq  [Jf  +  Jp+Jfr  j  +£o/io-^- 


dP 


dE 


VxB  =  fi o  [  Jf+  —  +  VxM  I  +£q11q— 


B  -  \  £o  E  +  P 

Vx  | - M  =Jf  + 

Mo 


dt 


.  B  ^  \  -  dD 

Vx  | - M  =  7/  +  — 

Mo 


dt 


(2.21) 


(2.22) 


(2.23) 


(2.24) 


(2.25) 
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B 

The  quantity - M  is  called  the  magnetic  field  and  denoted  H .  The  Ampere-Maxwell 

Mo 

law  for  macroscopic  media  is  thus 


Vx  #  =  //  +  —  (2.26) 

Replacing  (2.1)  with  (2.21),  and  (2.4)  with  (2.26),  we  arrive  at  Maxwell’s  equations  in 
matter 


V  D  =  pf 


(2.27) 


Vx£ 


SB 

~dt 


(2.28) 


VB  =  0 


(2.29) 


VxH 


(2.30) 


2.3  The  Wave  Equations 

2.3.1  The  Wave  Equations  in  Vacuum 

Let  us  consider  the  Maxwell  equations  in  vacuum  in  the  absence  of  sources.  Setting  p  —  0 
and  J  —  0, 


V£  =  0 


(2.31) 
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dB 

~dt 


(2.32) 


V  xE  —  — 


VB  =  0 


(2.33) 


VxB 


dE 


(2.34) 


Taking  the  curl  of  (2.32), 


VxVxl 


(2.35) 


Employing  the  vector  identity  VxVxf 


-v2y  +  v  (v-v 


VlE  +  V  V-E  = 


d  (VxB 


dt 


(2.36) 


Inserting  (2.31)  into  the  above  relationship, 


V2E 


dfvxBj 

dt 


(2.37) 


Inserting  the  sourceless  Ampere-Maxwell  law  (2.34)  in  vacuum  into  the  above  relationship, 


V2E  = 


(2.38) 


d2E 

~W 


~^—V2E 

£oMo 


(2.39) 
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Taking  the  curl  of  (2.34), 


V  x  V  x  B 


V 


dE 

x  £°^°~d7 


(2.40) 


-  V25  +  V 


dfvxE^j 

dt 


(2.41) 


Inserting  Faraday’s  law  (2.32)  into  the  above  relationship,  and  noting  from  (2.33)  that  V  ■ 

5  =  0, 


-V25 


(2.42) 


d2B 


-1— v2b 

£q  Mo 


(2.43) 


Recognizing  (2.39)  and  (2.43)  as  wave  equations  for  electric  and  magnetic  fields  which 

each  propagate  with  velocity  ——=  =  c,  which  is  experimentally  equal  to  the  speed  of 

V^oMo 

light  in  vacuum,  we  arrive  at 


d2E 

~W 


c2V2E 


(2.44) 


d2B 

~w 


c2V2B 


(2.45) 


2.3.2  The  Wave  Equations  in  Matter 

Taking  E  and  H  as  our  primary  fields,  Gauss’  law  in  matter  (2.27)  becomes 
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V  ■  (  £o E  +  Pj  —  pf 


(2.46) 


y  — _ — V .  p 

£o  £o 


(2.47) 


From  Faraday’s  law  in  matter  (2.28), 


V  x  E  =  —  po  - 


d[H  +  M 


dt 


(2.48) 


Equation  (2.29)  gives 


jUqV-  (H  +  M)  -0 


(2.49) 


V  H  =  —  V  M 


(2.50) 


The  Ampere-Maxwell  law  in  matter  (2.30)  may  be  expressed 


Vx  H  =  Jf  + 


d  (  £qE  +  P 


dt 


(2.51) 


As  was  the  case  in  vacuum,  we  assume  our  region  of  interest  is  free  of  sources  -  in  this  case, 
the  free  charge  density  and  free  current  density  pf  and  Jf,  respectively.  Taking  Pf  —  0  and 
Jf  =  0,  Maxwell’s  equations  now  reduce  to 


WE = -—VP 
£o 


(2.52) 
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(2.53) 


Vx£  =  —  jUq- 


d[H  +  M 


VH  =  —  V  -M 


(2.54) 


V  xH  — 


d  (  £qE  +  P 


(2.55) 


Once  again  similar  to  the  case  in  vacuum,  we  take  the  curl  of  the  Maxwell  equations  in¬ 
volving  the  curl  of  the  primary  fields.  Taking  the  curl  of  (2.48), 


VxVx£  =  —  jUqV  x 


d[H+M 


(2.56) 


Again  employing  the  vector  identity  V  x  V  x  V  =  —  V2V  +  V  •  V j ,  and  exploiting  the 

independence  of  the  temporal  and  spatial  partial  derivatives, 


-  V~E  +  V  ( V  •  E )  =  ~Ho 


dlVx  ( H  +  M 


(2.57) 


Inserting  (2.52)  and  (2.55)  into  the  above  relationship, 


2  ->  1  /  -A  d  d[£QE+P\ 

-  V-E - V(V-P)=-Uq—  — — - -  +  VxM 

e0  V  /  dt  dt 


(2.58) 


2  1  „  d2E  d2P  f  dM 

v^  +  -V  (V.P)  =  +Mo  V  x  — 


(2.59) 
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Forming  the  wave  operator  on  the  left,  and  identifying  c  = 


V£o  A*o 


as  previously, 


2  -*  1  d2E  d2P  /  dM\  1  /  - 

V2E~—  —  VX  —  VIV-P 


dt 2  ™  dt2 


dt  £q  V 


(2.60) 


Taking  the  curl  of  (2.55), 


VxVx//  =  Vx 


d  (  £qE  +  P 


dt 


(2.61) 


-V-H  +  y  (V-H)  =£o 


d  (  V  x  E  )  d  (  V  x  P 
+ 


dt 


dt 


(2.62) 


Inserting  (2.53)  and  (2.54)  into  the  above  relationship, 


V  H  —  V  ( V  -M )  =  £0^- 

dt 


-juo- 


dt 


+ 


d  [V  xP 


dt 


(2.63) 


d2// 


d2M 


Vz//  +  V(  V-Af)  =£0/io^  +  £oAto^^-Vx 


(2.64) 


Once  more  taking  c  =  and  forming  the  wave  operator  on  the  left, 

y£ofk) 


v2h 


1  d2H  1  d2M  „  dP  „ ,  „ 
-TT  =  -yTT-VxT-V(V'M 
c2  2  c2  dt2  dt 


(2.65) 


We  have  now  obtained  the  wave  equations  in  matter  for  £  and  H . 


2  1  d2E  d2P  („  dM\  1  |T7 

V-E-—^r^  =n  o^y  +  Ho  Vx  — - V  V-P 


;2  dt2  dt2 


df  /  £0 


(2.66) 
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(2.67) 


V2//- 


1  d2H 

c2  dt2 


1  d2M  dP 

-Vx- - V  V-M 


dt2 


dt 


This  set  of  coupled  hyperbolic  inhomogeneous  partial  differential  equations  will  prove 
resistant  to  the  common  methods  of  solution.  Following  an  overview  of  transformation 
optics,  which  may  be  viewed  as  a  method  for  circumventing  a  direct  solution  of  the  wave 
equations  above  for  a  desired  propagation  path,  we  will  return  to  these  expressions  in  an 
attempt  to  discover  an  iterative  solution. 
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CHAPTER  3: 

An  Overview  of  Transformation  Optics 


We  now  proceed  to  a  description  of  transformation  optics,  the  technique  currently  used  in 
the  design  of  metamaterial  structures  capable  of  redirecting  electromagnetic  radiation  as 
desired  through  those  structures. 

As  introduced  by  Pendry,  Schurig,  and  Smith  [2],  transformation  optics  employs  the  in¬ 
variance  of  the  Maxwell  equations  under  coordinate  transformations  to  convert  the  free- 
space  wave  solutions  in  a  coordinate-transformed  system,  referred  to  as  "electromagnetic 
space"  [9],  to  equally  valid  solutions  in  physical  space  with  electric  permittivity  and  mag¬ 
netic  permeability  as  determined  solely  by  the  coordinate  transformation  under  considera¬ 
tion.  A  full  understanding  of  the  elegance,  as  well  as  the  limitations,  of  the  transformation 
optics  approach  will  require  a  short  overview  of  coordinate  transformations  as  well  as  the 
covariant  formulation  of  classical  electrodynamics. 


3.1  Coordinate  Transformations 

3.1.1  The  Einstein  Summation  Convention 

Throughout  this  chapter,  repeated  indices  indicate  summation  over  all  components  up  to 
the  dimension  of  the  space  in  which  we  are  working.  That  is, 


n 

albi  =  £  aibl  (3.1) 

1=1 

where  n  is  the  dimension  of  the  space,  usually  3  (for  classical  mechanics)  or  4  (for  special 
relativistic  applications).  Furthermore,  in  order  to  avoid  confusion  between  contravariant 
tensor  components  and  exponents,  should  the  need  arise  to  express  a  tensor  component 
raised  to  a  power,  the  exponent  will  appear  outside  of  a  set  of  parentheses. 
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3.1.2  The  Jacobian 


Consider  an  invertible  coordinate  transformation  in  four  dimensions  from  some  coordinate 
system  xa  into  another  coordinate  system  xa'  given  by 


x  —x 


=  xv  (x°,xl,x2,x3) 


(3.2) 


and  similar  for  x0' ,  x2' .  x3' .  Coordinate  differentials  in  the  new  coordinate  system  are  given 
by 


,  dr«' 
dxa  =  -—dxa 
dxa 


(3.3) 


OL 

The  expressions  form  the  elements  of  the  Jacobian  matrix 


“  “  dxa 


(3.4) 


dxa  —Aaadxa 


(3.5) 


A  contravariant  vector,  with  components  denoted  Va  obeys  the  above  transformation  law. 


Va'  =Aa'aVa  (3.6) 

Just  as  the  coordinate  differentials  comprise  the  primordial  contravariant  vector,  the  gradi¬ 
ent  serves  as  the  prototypical  covariant  vector.  Note  that  for  some  function  f, 


3  f=AL  =  lL^  =  A<X  3  f 

a'f  dxa'  dxadxa'  a'  af 


(3.7) 


A  covariant  vector,  with  components  denoted  Va  will  obey  the  above  transformation  law. 
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3.1.3  The  Metric  Tensor 


The  covariant  metric  tensor  gjj  is  the  collection  of  inner  products  of  basis  vectors  in  a  given 
coordinate  system.  That  is, 


gtj  =  e/  •  e,-  (3.8) 

Note  that  due  to  the  commutative  nature  of  the  inner  product,  gjj  is  symmetric.  In  addition, 

ds2  —  gjjd.ddx'  (3.9) 

A  metric  gij  therefore  defines  a  geometry  on  the  coordinate  system  in  which  we  are  inter¬ 
ested  [16].  Covariant  vector  elements  may  be  related  to  their  corresponding  contravariant 
vector  elements  through  the  metric:  For  some  vector  V, 


Vi  =  gijVj  (3.10) 

Vi  =  gijVj  (3.11) 

These  operations  are  commonly  referred  to  as  “lowering”  or  “raising”  an  index,  respec¬ 
tively.  We  may  apply  these  relations  to  extract  a  useful  relationship  between  the  covariant 
and  contravariant  metric. 


Vi  =  gijVj  =  gijgjkVk  (3.12) 

5iVk  =  gijgjkVk  (3.13) 


Where  8k  is  the  Kroneker  delta. 
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(gijgJk-5hvt  =  0 


(3.14) 


Since  V  is  arbitrary, 


8ij8 


jk 


(3.15) 


and  we  conclude  that  gaP  is  inverse  to  gai 3 . 

3.1.4  Tensors 

A  contravariant  tensor  of  rank  2  is  a  geometric  object,  denoted  TaP  which  transforms  under 
the  coordinate  transformation  given  by  the  Jacobian  Aa'a  in  the  following  manner 

Ta'P'  =Aa'aA^TaP  (3.16) 

while  a  covariant  tensor  of  rank  2  Tap  transforms  as 

Ta'P'  =  AaaiA^p,Tap  (3.17) 

This  definition  may  be  extended  to  contravariant  and  covariant  tensors  of  higher  rank  by 
repeated  applications  of  the  Jacobian.  For  a  contravariant  tensor  of  rank  n, 


—  Aa  1  Aa~  Aa«  rra\a2---an 
Ci  1  ^2  •  •  • 


(3.18) 


while  for  a  covariant  tensor  of  rank  n, 


T  j  j  =Aa\Aa2,...Aan,Tc 

ala2--.an  a ^  n  < 


Ct\Cl2”'Cln 


(3.19) 


A  contravariant  pseudotensor  c£a^  of  weight  w  and  rank  n  is  an  object  with  components 
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which  transform  as  [16] 


=  |  A  | 


w &a\'  \al!  \an 
^  al"  a2 a 


TalQ 


2...dn 


(3.20) 


where  \A\  is  the  determinant  of  the  Jacobian.  Pseudotensors  of  weight  ±1  are  referred  to 
as  “tensor  densities.” 


3.2  Covariant  Electrodynamics 

3.2.1  The  Field  Tensors 

The  covariant  electromagnetic  field  tensor  in  Minkowski  spacetime  may  be  written 


/ 


Fap  — 


Q  fk  r2 


-=Z  B 


~BZ 

0 

Br 


By 


-Bx 

0 


(3.21) 


in  the  basis  (ct,x,y,z)  with  Minkowski  metric 


Sap 


(\  0  0  0  \ 

0-100 
0  0-10 
\0  0  0  -1  j 


(3.22) 


The  contravariant  electromagnetic  field  tensor  may  be  found  through 


/  0 


Fcci3=g«ngvlSF 


_& 

c 

0 


=*  B- 


%  -B, 


C 

~BZ 

0 

Bx 


_Ez\ 


B 


(3.23) 
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The  contravariant  magnetization-polarization  tensor  density  (weight  -1)  DJlaP  is 


(  o  Pxc 

PyC 

Pzc\ 

mccp  = 

O 

1 

—M- 

My 

-PyC  M, 

0 

—Mx 

\-PZC  - My 

Mx 

0  ) 

The  contravariant  displacement  tensor  density  (weight  - 

1)2>“0  [ 

o 

1 

£ 

—  DyC 

~Dzc\ 

X)aP  = 

Dxc  0 

~HZ 

Hy 

DyC  H- 

0 

-Hx 

yDZC  ~  Hy 

Hx 

0  / 

(3.24) 


(3.25) 


The  electromagnetic  field  tensor,  magnetization-polarization  tensor,  and  displacement  ten¬ 
sor  are  related  through 


3 


_L  paP  _  yj[aP 

Mo 


(3.26) 


3.2.2  The  Maxwell  Equations 

The  free  current  density  (weight  -1)  four- vector  is 

J}  =  (cp/,7/)  (3.27) 

From  Gauss’  Law  (2.27)  in  matter, 

V  D  =  pf  (3.28) 
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d\Dx  +  d2  Dy  +  <?3  Dz  —  Pf 


(3.29) 


^i-D01+^2-S02  +  <?3 -D03  =  -/?  (3.30) 

c  c  c  c  J 


d^01 +d2D02  +  d3T)03  =J°f 


(3.31) 


Since  £)00  =  0, 


doD00  +  <9iD01  +  d2®02  +  d3D03  -  Jf 


(3.32) 


While  from  the  Ampere-Maxwell  Law  (2.30)  in  matter, 


VxH 


(3.33) 


The  x-component  of  the  above  relation  is 


d2Hz  d3Hy  —  Jf,x  A  C(9q Dx 


(3.34) 


d2H-  -  d3Hy  =  J)  -  <90D01  (3.35) 

d2®21+d3®31  =j}-d0®01  (3.36) 

But  T)11  =0,  and  so 

di£u+<92®21+d3®31  =Jf-do®01  (3.37) 
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do^  +  dil^  +  ^D^+^S)31  =J} 


(3.38) 


Similarly,  the  y-component  of  the  Ampere-Maxwell  Law  in  matter  may  be  expressed 


d0D02  +  dlT>l2  +  d2T>22  +  d3  D32=J2f  (3.39) 

While  from  the  z-component, 


d0®03  +  diD13  +  <92D23  +  <932)33  -  Jj  (3.40) 

(3.32),  (3.38),  (3.39),  and  (3.40)  may  be  summarized,  in  the  Einstein  summation  conven¬ 
tion, 


d^v=J)  (3.41) 

We  note  that  the  above  expression  contains  both  Gauss’  law  and  the  Ampere-Maxwell  law 
in  linear  media.  Turning  to  the  two  sourceless  Maxwell  equations, 


VB  =  0 


V  xE  + 


dB 

~di 


0 


Note  that  (3.42)  may  be  expressed 


(3.42) 

(3.43) 


d\B\  -\-  <?2^2  T  <?3^3  —  0 


(3.44) 


d\  F23  +  <?2  F31  +  dj,Fi2  =  0 


(3.45) 
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While  due  to  the  antisymmetric  nature  of  the  electromagnetic  field  tensor  (3.21), 


d\F 32  +  d^F  13  +  c?3  F 21  —  0 


(3.46) 


From  the  x-component  of  Faraday’s  law, 


V 


x  E  + 


dB 

~di 


0 


(3.47) 


d^F^  —  <^3^2  +  cdoBi  —  0 


(3.48) 


cd2Fo3  +  C(hF 20  +  C  SqF 32  =  0 


(3.49) 


<hFo3  +  <hF  20  +  d{)F 32  =  0 


(3.50) 


Once  again  we  exploit  the  antisymmetric  nature  of  the  electromagnetic  field  tensor  to  obtain 


djF 30  +  <?3  Fq2  +  d{)F 23  —  0 


(3.51) 


Similarly,  from  the  y-component  of  Faraday’s  law  we  obtain 


d\ F 30  +  djFoi  +  do  F 13  —  0 


(3.52) 


d\  F os  +  d$F  jo  +  do  F 31  =  0 


(3.53) 


And  from  the  z-component  of  Faraday’s  law 
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d\  F02  +  <?2  F 10  +  do  F21  =  0 


(3.54) 


d\ F20  +  <?2  F01  +  d0F12  =  0  (3.55) 

Equations  (3.45),  (3.46),  and  (3.50)-(3.55)  contain  all  possible  three-element  cyclic  per¬ 
mutations  of  the  four  spacetime  indices,  and  may  be  encapsulated  as 

daFpa  +  dpF aa  4-  daF ap  —  0  (3.56) 


Or  equivalently, 


da 

:  (\e“»TSFrS 

)=0 

(3.57) 

where  eaPYs  is  the  four  dimensional  Levi-Civita  symbol.  The  quantity  -£a^5Fy§  = 

is  denoted  as  the  dual  tensor  to  the  electromagnetic  field  tensor. 

da^ap  = 

0 

(3.58) 

We  note  that  the  dual  electromagnetic  tensor 

(  0  Bx 

By 

Bz  \ 

jr®/ 1  = 

~BX  0 

—B  — 

Dy  c 

_E, 

c 

0 

Ey 

c 

_E* 

c 

(3.59) 

1 

Co 

1 

£* 

C 

0 ) 

1-  1 


may  be  obtained  from  the  components  of  FaP  through  B  — >■  E,  zf— »  =  B.  Equation  (3.58) 
contains  both  Faraday’s  law  and  V  ■  B  —  0.  Maxwell’s  equations  in  matter  may  thus  be 


28 


expressed  in  covariant  form  as 


d^v  =  Jj  (3.60) 

da^ali  =  0  (3.61) 


3.2.3  The  Covariant  Constitutive  Relations 

The  constitutive  relations  linking  the  primary  fields  to  the  material  response  fields  are  taken 
as 


P  =  £oXeE  (3.62) 

M  =  XmH  (3.63) 

for  homogenous,  isotropic,  linear  media,  where  the  scalar  values  Xe  and  Xm  are  the  electric 
and  magnetic  susceptibilities  of  the  material  under  consideration.  For  anisotropic  linear 
media,  we  take  these  susceptibilities  as  second  rank  tensors,  leading  to  the  relations 

P  =  e0xe*E  (3.64) 

M  =  Xm*H  (3.65) 

where  *  indicates  tensor  contraction.  By  components, 

P'  =  BoXeEj  (3.66) 
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(3.67) 


Ml  =  tlKj 

The  covariant  representation  of  these  relations  may  be  given  as  [4],  [10] 

m^v  =  X^vapFap  (3.68) 

where  are  the  components  of  the  rank  four  contravariant  susceptibility.  We  may 

conclude  from  the  known  nature  of  59T  as  a  tensor  density  of  weight  -1  and  the  tensor 
quotient  theorem  [16],  that  %  must  also  transform  as  a  tensor  density  of  weight  -1.  The 
traditional  permittivity  and  permeability  of  the  material  in  question  may  be  extracted  from 
the  susceptibility  tensor  density:  Inserting  (3.26)  into  (3.68)  [10], 

=  —  FaP  -  map  =  — Fa/3  -  (3.69) 

Mo  Mo 

Since  Fa$  is  antisymmetric,  it  may  be  expressed  [10] 

F “0  =  I  (F*  -  F?' “)  (3.70) 

^  =  j(*“Vv-/Vv)^v  (3.71) 

Inserting  this  relationship  into  (3.69), 

=  T;  (*“Vv- -  x“l>>‘v)  F„v  (3.72) 

Defining  Xa^8  as  [10] 
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(3.73) 


xccPy8  =  _L  fgaYgPS  -  g^rga5  - 

2  Ho  V 

Note  that  the  displacement  tensor  density  elements  are,  in  terms  of  this  new  susceptibility 
tensor, 


5)a^X“^vfjiv  (3.74) 

We  will  take  the  above  relationship  as  our  covariant  constitutive  relationship.  Note  again 
from  the  quotient  theorem  that  since  2)  transforms  as  a  tensor  density  of  weight  - 1  that  X 
must  also  transform  as  a  tensor  density  of  weight  -1. 

We  pause  for  a  moment  to  observe  that  our  constitutive  relationship  above  is  completely 
linear  and  therefore  does  not  explicitly  account  for  obserx’ed  nonlinear  effects  e.g.,  self- 
focusing,  frequency  doubling. 

We  may  extract  the  traditional  bianisotropic  coupling  coefficients,  such  as  the  permittivity 
and  permeability  from  the  elements  of  X.  From,  for  example,  the  (a/3)  =  (01)  term  [10], 

D°i  =X01^vFgV  (3.75) 

-cDx  =  Xm»vFgV  (3.76) 


-  cDx  =  XolooFoo  +  XomF01  +X0102Fo2  +  XomF03 

+  X01  WF10  +  XoinFn  +  X0U2F12  +  X0U3F13 

I  y'0120 jo  i  y'0121  r  i  v0122p  ,  y0123  r 

-\-A  t2oFA  C22FA  t23 

+  X0130F30+X013lF31  +  X0l32F32  +  X0133F33  (3.77) 
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_cD  =  ^x0102— -Xom—+Xolw— ~X0luB  +X0luB 

c  c  c  c  z  y 

+X0l20^y+X0l21B_x0l23Bx  +  x0l30^_X0l3lBy  +  X0l32Bx  (3.78) 

c  c 


cDx  =  (X0101  -X0110)  ^  +  (X0102  -X0120)  ^  +  (X0103  -X0130)  ^ 

+  (x0123  —  X0132)  Bx  +  (X0131  — X0113)  By  +  (X0112  —X0121)  Bz  (3.79) 


Comparing  this  relationship  to  the  common  bianisotropic  constitutive  relationship  [15] 


D  =  eE  +  ^H  (3.80) 

B  =  /iH  +  ^E  (3.81) 

where  £,  fi,  C,  are  the  spatial  tensors  which  couple  the  electric  displacement  and  mag¬ 
netic  fields  to  the  electric  and  magnetic  fields, 

Elj  =  (X01°2 -X0lj0)  (3.82) 

Note  that  since  the  electromagnetic  field  tensor  F  is  antisymmetric,  our  constitutive  rela¬ 

tionship  (3.74)  requires  that  Xa^vbe  antisymmetric  with  respect  to  the  pairs  /iv.  There¬ 
fore 


£lj  =  4 X0l0j  (3.83) 

Similar  calculations  yield  expressions  for  the  remaining  components  of  the  permittivity 

[10] 
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(3.84) 


£lj  = 


3.3  Transformation  Optics 

Having  constructed  the  contravariant  rank  4  modified  susceptibility  tensor  density  Xa^5 
[15],  and  analyzed  the  transformation  properties  of  tensors  [16],  we  may  now  investigate 
the  transformation  properties  of  the  material  parameters  contained  within  the  susceptibility 
tensor.  As  a  contravariant  tensor  density  of  rank  4  and  weight  -1,  under  a  coordinate  trans¬ 
formation  described  by  the  Jacobian  A  aa ,  the  contravariant  components  of  the  redefined 
susceptibility  X  transform  as  [10] 

Xa'P'y8'  =  |  A\~l  AaaAPpAy,7As'sXaprS  (3.85) 

We  now  investigate  the  way  in  which  the  Maxwell  equations  in  media  transform  under  such 
a  coordinate  transformation.  From  (3.60)  and  (3.61), 

dfj.X )tlv=Jj  (3.86) 

da^ap  =  0  (3.87) 

Applying  the  coordinate  transformation  given  by  the  Jacobian  with  elements  Aa'a  to  the  left 
hand  side  of  the  Gauss-Ampere-Maxwell  law  (3.84), 


=  IAP1  A"  =  HP 1  A^A^Jf 


(3.88) 


But 
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,  An'  =  W_M?_  =  W  = 
V-'  ,U  rXi'A  r)x^  dx^ 


(3.89) 


and  so 


d^V  =\A\~lAVyj} 


(3.90) 


Since  the  four-current  density  is  a  tensor  density  of  rank  1,  weight  -1,  it  will  transforms  as 


Jf  =  I 


(3.91) 


The  Gauss-Ampere-Maxwell  law  therefore  transforms  as 


d^v  =JJ 


(3.92) 


and  is  thus  form  invariant  under  coordinate  transformations.  Turning  to  the  Faraday-Gauss 
law, 


=Aaa,Aa'aA^da^  =0 


(3.93) 


All  of  the  Maxwell  equations  are  therefore  form  invariant  under  coordinate  transforma¬ 
tions,  as  we  might  expect  of  such  a  set  of  tensor  relationships.  A  valid  solution  set 
,  T)“^,  Xa^5  to  Maxwell’s  equations  in  one  coordinate  system  will  provide  a  valid 
solution  set  ,  T)a ,  Xa'^'^8'  in  another  coordinate  system  according  to  our  transfor¬ 
mation  properties  [10] 


=Aa'aA^aP 


(3.94) 
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T)a'P'  =  \A\~X  AaaA^pT)aP 


(3.95) 


xa'pys'  =  iA\-1AafaAp'pA7'yAs'sXa^s  (3.96) 

The  form  invariance  of  the  Maxwell  equations  and  the  resulting  transformation  properties 
of  the  modified  susceptibility  tensor  density  X  is  the  concept  at  the  heart  of  Transformation 
Optics. 


3.3.1  Spatial  Transformations 

Turning  now  to  the  special  case  of  transformations  which  do  not  involve  the  time: 

t’  =  t  x'=x\x,y,z)  y'  =  y\x,y,z)  z'  =  z\x,y,z)  (3.97) 

Note  that  under  this  set  of  restrictions 

A°q  =  1  A0-  -  0  (3.98) 

for  all  i  from  1  to  3  (spatial  index).  From  (3.84)  the  components  of  the  electrical  permittivity 
transform  as 


//  =  -|x°W/  =  |A|-1A0^A0^/  xapyS 


(3.99) 


From  (3.98),  the  only  a  or  y  which  survive  the  summation  are  a  —  0  and  y  =  0. 


//  =  A  lAr^^A^X0^05 


(3.100) 


But  from  (3.84),  —X°PQS  =  e^5 .  Reindexing,  [10] 
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(3.101) 


e'V  =  \A\~l  A1' :Aj' j£ij 

3.3.2  Application:  Cylindrical  Cloak 

As  our  first  investigation  into  the  design  of  an  electromagnetic  redirection  structure,  we  fol¬ 
low  the  development  of  a  right  circular  cylindrical  device  in  [9].  Let  unprimed  coordinates 
represent  physical  space  and  primed  coordinates  represent  our  electromagnetic  space.  We 
consider  a  coordinate  transformation  of  the  form 

r=  — - r'  +  a  (3.102) 

b 

</>  =  </>'  (3.103) 

z  =  z'  (3.104) 

ct  —  ct'  (3.105) 

in  cylindrical  coordinates,  where  a  and  b  are  positive  real  numbers,  with  b  >  a.  Such  a 
transform  maps  the  line  r'  =  0  onto  the  circular  cylinder  r  =  a  and  the  circular  cylinder  r'  = 
b  onto  the  circular  cylinder  r  —  b.  Figure  3.1  illustrates  such  a  coordinate  transformation. 

Let  us  assume  vacuum  in  the  electromagnetic  space.  Then  £l'J'  —  £q 81-,  and  fl1' i'  = 

As  the  transformation  (3.102)  does  not  involve  the  time,  we  may  apply  (3.101)  in  order  to 
determine  the  transformed  permittivity  in  physical  space 

£ij  =  \A\-lAii,Ajr£i'f  (3.106) 

The  Elements  Aa  ,  of  the  Jacobian  are 
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Figure  3.1:  Transformation  for  Cylindrical  Cloaking  Structure.  Source:  [9] 


,«  _  ^ 


(3.107) 


And  so 


4° 

A  0, 


d  ( ct ) 

djrf) 


(3.108) 


d(ct) 

dx1' 


-0 


(3.109) 


A' 


o' 


dx1 

djrt) 


0 


(3.110) 


for  all  spatial  indices  z,z-/  =  1,2,3-  Computing  the  remaining  elements  of  this  Jacobian, 


1  dx  d(rcos(<j>)) 

v  dx'  dx' 


(3.111) 
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(3.112) 


A1 


l' 


d  (rcos(0))  dr'  d  (rcos  (0))  d(/)'  d  (rcos  (0))  dz! 

dr'  dx'  d(j)'  dx!  dz'  dx' 


A1 


l' 


d_ 

dr' 


b  —  a 
b 


r'c os  (0;)  +c/cos  (0') 


dr'  d  (rcos  (0'))  d(f>' 
dx'  d(j)'  dx' 


Let  7? 


&  —  a 
b 


Alv  =  R cos  (0/)  —  rsin  (0;)  ^ 

O'  A 


<9jc' 


<9  ( v  y2 + y2 

A!j,  =  7? cos  (<^/)  — 2 — - —  —  rsin  (</) 


d  [  tan  1  ( ^ 
dx' 


Alv  =  7? cos  (0')  ^  +  rsin  (0') 


*/2(l  +  ^ 


A 1 1,  —  Rcos2  (</)  +  —f  sin2  (</) 


l  _  d*1  _  d  (rcos  (0)) 
“  (9a'2/  dy' 


A 


l 

2' 


a  r  , 

^—7  \Rr  cos 
dy  L 


(0;)  +acos  (</)] 


Al2,  —  (Rr'  +  a) 


d  (cos  (0')) 
dy' 


+  7? cos  (0;) 


d^_ 

dy' 


(3.113) 

(3.114) 

(3.115) 

(3.116) 

(3.117) 

(3.118) 

(3.119) 

(3.120) 
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A 


l 

2' 


(Rr  +  a) 


+  Rc os  (07) 


J  r)  ( \Jxa  +/2)  d  (  a Jx'2+y'2 

J4»2  =  -(^  +  a)2_  V  ^+flcos(f)  ^ 


.vV 


AV  =  -(*^0)^+*co.(^ 


a'2, = -  (s/+a)  <^mpn  +Rcos  m  sin  (f ) 


A  '2.  =  —  ^  cos  (^')  sin  ((/>')  =  f  cos  (0')  sin  (0') 


A1 2 :  =  —  —f  j  cos  (0;)  sin  (0;) 


A 


l 

3' 


(9.x:1 

a^7 


a* 

a? 


=  o 


A\'  =  AfAlAA  =  A  [(Rr’+a)  sin  (0')] 
A\=(Rr'+a)3-^pi+Rsi n(0')% 


( Rr '  +  a ) 


dxr 


+  R  sin  ((j)') 


dr1 

dx! 


(3.121) 

(3.122) 

(3.123) 

(3.124) 

(3.125) 

(3.126) 

(3.127) 

(3.128) 

(3.129) 

(3.130) 
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A2i,  -  -  (Rr'  +  a)  l^^+Rs in(0/)|^  (3.131) 

A2,  =  (- 

A2i'  =  -7  sin  (0')  ^7  =  sin  (3.1.33) 

A2t,  =  -  r -  cos  ((/)  sin  (<)/)  (3.134) 


.  ,  x  sin  (</>') 

(i?r  +  a) - )-A-  4-  7?  sin 


(3.132) 


A2v  —  (ft—  -7)  cos  (</)  sin  (</)  (3.135) 

a2?  =  a(y»  =  9_  [{R,.’  +  a)  sin  (f )]  (3.136) 

A^  =  (^  +  «)^A^si»(f)|  (3.137) 

dtt)  d/ 

A22'  =  (RtJ  +  a)  — ^y—  +  R sin  ((/)')  —  (3.138) 

A22'  =  (R,J  +  a)  (y~p2gy  +  ~l^j  +  Rsin(^)gy  (3.139) 

A22,  =  (1 Rr'+a )  +  ^  +i?sin (^)/) (3.140) 
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^,=  («/  +  a)1-Sin,2(f)+/es  in2(f) 

(3.141) 

AJ2,=  (&'  +  a)COs2,W,)+Ssi„J(f) 

(3.142) 

A22i  —  Rsin2  ((j)')  +  —j cos”  [(j)') 

(3.143) 

(3.144) 

The  remaining  components  of  the  Jacobian  are  somewhat  simpler.  For  i' 

=  1,  2,  we  obtain 

A2  -  3Z  -0 

A'~  Ac'-0 

(3.145) 

While 

A3  -  dz  -  1 

A  3'  "  dz'  ~  1 

(3.146) 

Our  Jacobian  therefore  takes  the  form 


/  Rcos2  ((j)')  +  p-sin2  ((j)')  (R—  y)  cos  (0;)  sin(0;)  0^ 

A—  (R—fr)  cos  ((j)')  sin  ((j)')  Rsin2  (0')  +  p-cos2  ((j)')  0 

\o  0  1/ 

The  determinant  \A  \  of  this  Jacobian  is 


(3.147) 
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\A\  —  (rc os2  (0;)  +  -^sin2  (0')  j  sin2  (0;)  +  -^cos2  (0;)  j  —  (^R  cos2  (<K)  sin2  (Q') 

(3.148) 


\A\  —  ^7?2  +  cos2  (0')  sin2  (<}>') +  ~y  (cos4  (0')  +sin4  (0'))  —  (r~  cos2  (0')  sin2  (^') 


(3.149) 


|A|  =  ^  (cos4  ((/)')  +sin4  ((/>'))  +2^- cos2  (</)  sin2  (0')  (3.150) 


lAl  =  ^(cos2(f)+sin2(0/))^ 


(3.151) 


lAl  =  ~t 

r 


(3.152) 


Applying  our  transformation  relationship  (3.106)  for  purely  spatial  transformations,  with 


£ij'  =  £o8ji, 


£lJ  =  \A\  lA\,AJf£ o8lr 


(3.153) 


£ij  =  \A\  1AijlAJf£ o 


(3.154) 


£ij  =  £o  \A\Alj,  (AT) j. 


(3.155) 


£ij  =  £o\A\~l  (AAt)’j 


(3.156) 
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(3.157) 


£  =  £o  | A  |  1  AAt 


£  =  £0Ri- 


(  Rcos2  (00  +  psin2  (00  (ft— p/c os  (00  sin  (00  0^ 

(ft  — p)  cos  (00  sin  (00  7? sin2  (0')  +  pcos2  (00  0 

0  0  1/ 

^ftcos2  (0;)  +  7  sin2  (00  —  7)  cos  (00  sin  (00 

(ft  -  p)  cos  (00  sin  (00  ft  sin2  (00  +  7  cos2  (00 
0  0 


x 


V 


ell=e°^:((^cos2(0/)  +  ^sin2 


+  R-^ 


cos" 


(00 si 


snr 


£l1 = £ok-  [r2cos4  (^0 + 72 sin4  (00  +  ( + 62 ) cos2  W) si 


r 

72 


sin 


.11 


=  £o~"  (  ft2  cos2  (00  (cos2  (00  +sin2  (00)  +  -72  sin2  (00  (sin2  (00  +cos 


£ 


11 


(3.158) 


(3.159) 


(3.160) 


(3.161) 

(3.162) 
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g12  =  £o^  (r  (R  ~  cos3  (0;)  sin  (f )  +  £  {R  ~  £)  cos  (0')  sin3  (0') 

+  i?  (*  “  ?)  C°S  W  Sin3  ^  +  £  (R  ~  ?)  CO$3  W  Sin  W  )  (3'163) 

e  12  =  £o ^(r(r-^)  cos  ((/)')  sin  ((/)')  +  £(*-£)  cos  (0;)  sin  (0') )  (cos2  (0;)  +  sin2  (0;) ) 

(3.164) 

el2  =  eo^(^+^)  (i?-^)cos(f)sin(f)  (3.165) 

£12  =  £o^  (V  -  cos  (0')  sin  (f)  (3.166) 

£13  =  0  (3.167) 

£21  =  £o^  (i?  -  cos3  (f)  sin  (0')  +  £(*-£)  COS  (0')  sin3  (0;)  + 

7?  ^/?  —  cos  ((/)')  sin3  ((/)')  +  —t  (/?  —  -^)  cos3  (0;)  sin  (0;)  ^  (3.168) 

£21  =  £o^;  (^-  ^7)  cos  (0;)  sin  (0;)  +  r—  [r-  cos  (0;)  sin  (0;)  j  (cos2  (0')  +sin2  (0')) 

(3.169) 

£21=e°^-(W^)  (/?-£)  cos  (0')  sin  (0')  (3.170) 
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£22  =  £q  y—  cos2  sin2  (00  +^2sin4  ((j)')  +  pr cos4  (</>')^J  (3.173) 


£22  =  £o^-  ^2sin2  [(j)')  (cos2  ((j)r)  +sin2  ((/>'))  +  -^-cos2  ((/)')  (sin2 

(0)+COS2(f))^ 

(3.174) 

e22 = £ok  (■ Rlsin' 2  (*') + ?2 cos2  (*')) 

(3.175) 

£31  _  £32  _  q 

(3.176) 

£33  =  1 

(3.177) 

Collecting  the  elements  of  £, 

^  /^2cos2(0')  +  ^sin2^')  (ft2  — cos  (</>')  sin  (</)')  0^ 

£  =  £o—  ^i?2  —  ^  j  cos  (0')  sin  (^')  tf2sin2  (0')  +  72  cos2  (0')  0  (3.178) 

V  0  0  1) 
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in  accordance  with  results  claimed  in  [9].  Note  that  since  r'  — »  0  as  r  — >  a.  the  susceptibil¬ 
ity  components  j}1  ,y}2 ,%21 ,  and  y22  are  expected  to  diverge.  This  divergence  indicates  a 
possible  issue  which  may  arise  when  performing  numerical  computations  in  order  to  deter¬ 
mine  the  response  fields  (see  Chapter  5).  Specifically,  we  anticipate  the  possibility  that  the 
value  attained  from  those  calculations  may  vary  considerably  depending  on  the  fineness  of 
our  mesh. 

Converting  to  cylindrical  coordinates  in  the  basis  (r,  <j>,z), 

£'■'/  =  |  A\~x  A1' ft' j£ij  (3.179) 

(V7)  =  \A\~l AEcanesian  ( A)T  (3.180) 

where  A  is  the  Jacobian  from  Cartesian  to  cylindrical  coordinates,  given  by  the  transforma¬ 
tion 


r  =  \J  x 2  +y-  ,  (j)  —  tan 


-i  fy 


z  —  z 


(3.181) 


Computing  the  Jacobian, 


djf_ 

dxJ 


—>A  = 


( ita. 

dx 

d(j) 

dx 

\  dx 


dr  dr  \ 
dy  dz 
dtp  dtp 
dy  dz 
dz  dz 
dy  dz ) 


(3.182) 


/  - 
r 

Z 

r 

°) 

/cos (0) 

sin(^>) 

o\ 

zZ 

x_ 

» 

—  sin  (0) 

cos  (<p) 

0 

r2 

r2 

r 

r 

\o 

0 

l) 

V  o 

0 

1/ 

(3.183) 


The  permittivity  in  cylindrical  coordinates  is  thus 
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i'j' 


/cos  ((f))  sin  (<f>)  0\ 

sin  (0)  cos  (0)  q 

r  r 

o  01/ 


^ R2  cos2  (0)  +  ^  sin2  (0) 

(r2- 

J 

I  cos  ((f))  sin  (</)) 

0^ 

(^2  —  pi)  cos  (0)  sin  (cj)) 

R2sin 

2  ((/>)  + COS2((/)) 

0 

V  o 

0 

1/ 

/cos (0) 

sin((/>) 

r 

(A 

x  sin  (<j>) 

COS  (0) 
r 

0 

(3.184) 

V  o 

0 

V 

jf  j' 

£  J  )  = 


r' 

£«r 


/R2c  os(0)  R2sm((j))  0\ 


i(0) 


r- 

0 


s(0) 


0 


0 

1/ 


sin(0) 

r 

COS  (0) 
r 

0 


o\ 

0 

1/ 


(3.185) 


£*'/  I  = 


r' 

£°r 


( R 2  0  0\ 
0^0 
\0  0  \) 


(r'R  0  0\ 


=  £o 


0  7R 


\0 


0 


o  i/ 


(3.186) 


0 


£o 


0  0 


°) 

0 


r— a  / 
R2  / 


(3.187) 


Note  that  these  results  are  not  in  accordance  with  [9],  although  both  results  are  diagonal  in 
cylindrical  coordinates.  We  attribute  the  discrepancy  (of  a  factor  r)  to  the  author’s  neglect 
of  the  tensor  density  nature  of  the  permittivity  and  subsequent  omission  of  the  necessary 
factor  \A  |  1  —  r.  Regardless,  the  permittivity  is  diagonal  in  cylindrical  coordinates,  as  we 
might  have  predicted  based  upon  symmetry  requirements.  Again  we  see  that  at  least  one 
component  of  the  susceptibility  is  expected  to  diverge  at  r  —  a. 
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CHAPTER  4: 

An  Iterative  Approach  to  the  Wave  Equations 


In  this  chapter  we  will  discuss  an  alternative  approach  to  the  wave  equations  in  matter 
as  determined  in  chapter  two.  We  will  then  consider  this  approach  when  applied  to  the 
problem  of  scattering  of  a  pure  plane  electromagnetic  plane  wave  from  a  structure  of  known 
geometry  and  electromagnetic  susceptibilities. 
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frequency  domain. 


V2E  +  k2E 


-—P-i(OU0  (v  xM]  -  —V  (v-P 

£0  V  )  £0  V 


(4.5) 


where  k  =  —  and  we  have  used  the  relationship  c  = - . 

c  v/eoMo 

frequency  Fourier  transform  of  (4.2), 


Similarly,  from  the  time— > 


(0 


CQ- 


V^H  +  ^rH  = TM  +  ia>  (  V  x  P 

c~  c~ 


V  V-Af 


(4.6) 


V2fi  +  k2fi  =  -k2lti  +  i(Q  (VxP)  -v(v-itf) 


(4.7) 


where  k  is  again  equal  to  — .  We  recognize  (3.5)  and  (3.7)  as  a  set  of  inhomogeneous, 

c 

coupled  Helmholtz  equations  for  the  components  of  E  and  H. 


4.2  Constitutive  Relations 

We  wish  now  to  introduce  a  set  of  constitutive  relations  between  the  polarization  and  mag¬ 
netization  fields  P  and  M.  and  the  electromagnetic  fields  E  and  H .  respectively.  In  linear, 
isotropic  media,  the  relations 


P  -  £oXeP 

M  =  %mH 

where  xe  and  X"‘  are  the  electric  and  magnetic  susceptibilities,  and  are  employed  to  de¬ 
scribe  the  response  of  the  polarization  and  magnetization  to  the  fields.  In  order  to  account 
for  the  possible  anisotropy  of  our  media,  we  express  each  of  the  susceptibilities  as  the 
second-rank  tensors  xe  and  X'n ■  The  constitutive  relations  now  read 
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(4.8) 


P  =  £oXe*E 

M  =  xm*H  (4.9) 

where  *  denotes  tensor  contraction.  By  components, 

f«=foE  <4'10> 

0-1 

<4-H) 
P  =  1 

Our  Fourier-transformed  wave  equations  now  read 

V2I  +  k2E  =  -k2  [xe  *E^j  -iajiM o(Vx  (*"'*#))  -V(V-  (V  *£))  (4.12) 

V2£  +  k2fr  =  -k2  (xm  *  #  )  +  i0)£0  (v  x  (V*I))-V(V-  (xm*fi))  (4.13) 

4.3  Incident  and  Response  Fields 

We  now  consider  separately  the  incident  fields  and  the  material  responses  to  those  fields. 
That  is, 

i  =  P0  +  ir  (4.14) 

B  =  fi0  +  fir  (4.15) 
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where  Eq,  Hq  are  the  incident  fields  and  Er,  H,  the  fields  due  to  the  material  response.  From 
(4.12), 

V2  +  £,-)  +k2  (4  +  ^-)  = 

-  k2  [xe  *  (4)  +  i)  )  -  ico/M)  (V  x  [xm  *  (#o  +  fir)  )  )  -  v  (v  •  (t  *  (4  + 1)  )  ) 

(4.16) 

Let  us  model  the  incident  fields  as  arriving  through  vacuum.  Then  the  incident  fields  are 
expected  to  obey  the  sourceless  Helmholtz  equations: 

V2Eq  +  k2Eo  =  0  (4.17) 

Thus, 

v2Er+k2Er  = 

-  k2  (V  *  (e0  +  ir)  )  -  icon 0  (v  X  [xm  *  (hQ  +  fir)  )  )  -  V  (v  •  (V  *  (4  +  ir)  )  ) 

(4.18) 

Similarly,  from  (4.13) 

V2fir  +  k2fir  = 

-  k2  (xm  *  (4  +  fir)  )  +  /we0  (v  x  [xe  *  (A  +  i)  )  )  -  V  (v  •  [xm  *  (4  +  fir)  )  ) 

(4.19) 

Note  the  near  symmetry  in  the  source  terms  of  the  above  two  equations  for  the  response 
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fields.  We  define 


Fi(x,F)  =  -k2(xe*X^-ia)^(vx  (z"'*?))-v(v-(V*x))  (4.20) 

f2  (x,  y  )  =  -k2  (xm  *  x)  +  irn 0  (vx(f*y))-v(v.  (*m  *  x'j )  (4.2 1 ) 

Note  that 

Fi  (xl+X2,Yi+Y2)  = 

-  k2  (xe  *  (xi  +X2)  )  -  im o  (v  X  (xm  *  (fi  +  f2)  )  )  -  V  (v  ■  (>  *  (f  ,  +  X2)  )  ) 

(4.22) 

=  - k 2  (>  *Xi)  -  k2  (V  *X2)  -  ico/M ,  (v  x  (Vw  * ?i) )  -  /como  (v  x  (x'n  *  *2)  ) 

-v(v.(^*X!))-v(v.(^*x2)) 

F,  =  F\  +Fi  (x2,?2)  (4.23) 

And  similar  for  F2.  From  (4.18)  and  (4.19), 

V2Fr  +  ]cEr  =  Fi  (f0, 4)  +  Ft  (Fn A)  (4.24) 

V2/?r  +  k2Hr  =  Fi  (S0  A)  +  F2  (ArX)  (4-25) 
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4.4  An  Integral  Equation 


The  Green  Function  for  the  Helmholtz  equation  in  two  dimensions  is  known  (see  Appendix 
A)  to  be  given  by 


Gk o  {r,r')  =  l~H^  ( k0  |r-r'|)  (4.26) 

For  the  Helmholtz  operator  V  +  Icq,  where  HQ  is  the  first  Hankel  function  of  the  first 
kind.  With  this  Green  function  in  hand,  we  may  express  the  solution  to  our  inhomogeneous 
Helmholtz  equations  (4.24)  and  (4.25)  as 


Er  =  Gk  (r,  7)  (V,  (e0  (7) ,  H0  (7) )  +  F\  (ir  (r1)  ,  Hr  (r;) )  )  d2r'  (4.27) 

Sr  =  f~  Gk  (r,  7)  (f2  (4  (7 )  Jo{7))+  F2  (fir  (/) ,  ir  (r) )  )  d27  (4.28) 

We  have  now  converted  our  differential  equations  for  the  response  fields  in  equations 
(4.24)  and  (4.25)  into  a  set  of  coupled  integral  equations  for  these  fields.  In  an  attempt  to 
keep  our  expressions  under  control,  we  further  introduce  the  notation 


F[(X,y)  =Fl(X(r),Y(r) 


(4.29) 


and  similar  for  F2.  Our  integral  equations  for  the  response  fields  now  read 


Er=  /  Gk{r,7)F f  (Eq,Hq)  d27  +  /  Gk  (r./)  F{  (EnHr)  d 


27 


(4.30) 


Hr  -  /  G,  (r, 7)  F{  (Hq,Eq)  d27  +  /  Gk  (r, r')  F|  (Hr,Er )  d 


27 


(4.31) 
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4.5  The  Neumann  Series 


We  propose  to  attain  an  approximate  solution  to  equations  (4.24)  and  (4.25)  for  the  response 
fields  through  the  following  iteration  method:  Insert  the  entirety  of  equation  (4.27)  and 
(4.28)  for  Er  and  Hr,  respectively,  back  into  the  last  terms  of  (4.27)  and  (4.28).  This  method 
produces  what  is  known  in  mathematics  as  a  Neumann  series,  which  is  often  employed  in 
the  approximate  solution  of  integral  equations  [17].  The  technique  is  also  commonly  used 
for  quantum  scattering  problems,  where  truncation  of  all  but  the  first  term  of  the  resulting 
series  is  known  as  the  Born  approximation  [18]. 


P-rJ 


Er  =  I  Gk(r./)F[  d 

+  rGk(r/)F?(  (  r  Gk(r',r")Ff  (e0,H0)  d2r"+  I  Gk(r',r")F f  (EnHr)  d 


Gk  (r1  ,r")  fT  (Ho,E0)  d2r"+  I  G,  (r',?")  FI"  (Hr,Er)  d 


cr" 


r2-*// 

v  r 


p-fj' 


d2r  (4.32) 


Er  =  f^Gk  (?,?')  F{'  (g0 ,4)  d2r' 


/oo 

Gk  r  )  Ff 

-oo 

/oo 

Gk  r  )  F[' 

-oo 


d27",f  Gk  (??')  F?  (h0,E o)  d2r'^j  d2r' 
(ErM,)  d2r'\j^  Gk(r’,r")F?  (ftr,£r)  d2r'^j  d2r' 


(4.33) 


where  we  have  used  the  linearity  of  the  function  F\  from  equation  (4.23)  in  the  last  step. 
Provided  the  convergence  of  the  resulting  series,  we  may  continue  this  iteration  process 
as  we  wish  in  order  to  attain  an  approximate  solution  for  Er.  Similarly,  for  the  response 
magnetic  field  Hr, 
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4  =  f  Gk  {7,7')  F{  (4,4)  d2r' 


i2->n  , 
a  r  + 


5 


G,  {7'. 7")  F(  (44) 


J2r"  + 


G/c 


r'r" 


)n 


F  FI 

*^ri  llr 


d2r'  (4.34) 


4  =  /  Gk  {7,7')  F{'  (fio,£0)d 

J  — oo 

/oo 

G*  (?/)*!' 

-oo 
poo 

+  Gk(r,r')FlJ 


;2p/ 


Gk  (r,  r")  4"  (tf04o)  dV,  /  G,  (/,?")  if  (44)  dV 


Gk  (r;,  r")  Ff  ( Hr,Er )  J2r",  /  Gfc  (r;,  r")  if'  ( Er,  Hr )  d2r"  J2/ 


(4.35) 


Identifying  the  first  term  of  (4.33)  and  (4.35)  as  the  first-order  response  in  %,  second  term 
as  the  second-order  response,  and  the  final  term  the  remaining  higher-order  response,  we 
may  continue  the  iteration  process  as  we  desire  by  inserting  each  of  (4.33)  and  (4.35)  into 
the  expressions  for  the  response  fields  of  each  of  these  expressions,  each  time  attaining  the 
next  higher  order  solution  to  our  wave  equation. 


4.6  Incident  Plane  Wave 

As  a  precursor  to  our  implementation  of  the  iteration  scheme  outlined  above  (see  Chapter 
5),  let  us  consider  the  material  response  of  a  two-dimensional  structure  embedded  in  the  x-y 
plane  when  subject  to  an  incident  electromagnetic  plane  wave,  polarized  in  the  z  direction, 
with  angular  frequency  (Oq,  travelling  in  the  x,  direction,  and  arriving  through  vacuum. 


F0  = 


(4.36) 
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where  ko  is  the  wave  vector,  oriented  in  the  direction  of  propagation,  with  magnitude  given 

u  ,  0)0 

by  k0  =  — . 

c 


4.6.1  First  Order  Response 

The  temporal  Fourier  transform  Eq  (7.  to)  with  respect  to  waves  propagating  outward  is 

,  poo 

E0=  Eoe^-^e1™  dtz  (4.37) 

J  — oo 

}  poo 

E0  =  E0eikoX  e^-^dtz  (4.38) 

J  —oo 


E0  =  2 7lE0eikoX8  (oj- (0o)z  (4.39) 

where  8  (ft)  —  (Do)  is  the  Dirac  delta  distribution.  For  plane  electromagnetic  waves,  Hq  is 
determined  by  our  choice  of  polarization  and  propagation  directions  to  be 


Hq  =  -Hoe^-^  y 


(4.40) 


Where  Hq 
wave  is 


Eq 

H-QC  Y  /i o 


EqcEq.  The  temporal  Fourier  transform  of  this  incident  plane 


Hq  =  -2nH0eikoX8  {(0-(0o)y  (4.41) 

We  now  take  %e  =  x'n  in  order  to  ensure  that  the  structure  is  impedance  matched  to  free 
space,  eliminating  reflection  -  a  desired  property  in  a  CDEW  application.  The  first-order 
response  electric  field  is  thus,  from  (4.33), 


Erl  = 


Gk  {r,r')F{'  (ioM  d 


(4.42) 
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ir-'=  !-„L4H«)(k 

x  {-k2  (x  *Eo  (r;) )  -  i(QH0  (' V  x  [x  *H0  (r;) )  )  -  V  ( V  ■  (x  *  E0  (/) )  )  )  d2r’  (4.43) 

x  ( k 2  (xe  *  4)  (r ) )  +  /COMO  (v  x  (*,„  *  4  {r ) )  )  +  V  (v  ■  (xe  *  4)  {r ) )  )  )  d2r'  (4.44) 

Inserting  the  expressions  (4.39)  and  (4.41)  for  Eq  and  Hq,  respectively,  and  exploiting  the 
relation  Xe  =  Xm  for  nonreflectance, 

4,1  =  -8(CD-CDo)™  J  H[q]  (k\r-r\) 

(k2  [x  *  (E0eikoX'  z)  )  -  i(Oii0  (vx(p  H{)eik{)K'  y)  )  +  V  (v  •  [x  *  ( z)  )  )  ) 

(4.45) 

4,1  =  -5(®-(no)y  J  (k\r-r'\) 

(k2EQ  (x  *  (eikoV  z)  )  -  iffl/xoffo  (Vx(p  y) )  +  £0 V  ( V ■  (* *  (V*0*'  z)  )  )  )  d  V 

(4.46) 


Let  us  express  %  in  the  basis  {x,y,z}  as 
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(X 11  X12  0 

X  =  X21  X22  0 

V  0  0  X33 


(4.47) 


since  all  of  the  TO-derived  susceptibilities  under  consideration  in  this  work  take  the  above 
form.  Furthermore,  since  all  of  the  components  Xij  have  no  z-dependence,  the  tensor  con¬ 
tractions  of  equation  (4.46)  are 


X*(eik°xz)=eik°xX33Z 


(4.48) 


V  ■  (x*eikoXz)  =  V  ■  (eik°xX33z)  = 


d  ( eik°xX33 ) 


(4.49) 


iko-x  A  _  Jk()X  ^X33  * 


V-  (  x*elK0Xz)  =  elK°x-^z 


(4.50) 


As  stated,  X33  does  not  vary  with  z,  and  so 


V-  (x*eikoXz)=  0 


(4.51) 


Turning  to  the  curl  term  in  (4.46), 


X  *  eik°xy  =  eik°x  {XvJ  +  X22 $) 


(4.52) 


[” ^7  (  (  iknxA  M  (elk°X  (xn$l ,1  +X2281 ,2)) 

v  x  (**  (ek°  yj  J J  =  eiW— - 4 - _1 - lAT 


(4.53) 
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Vx(p(Alf)) 


d  (eikQXX2i)  \  jg+  d  (e**Xn )  |  /  d  {eihxX22)  _  d  ( eik°xXn )  \  , 

<9z  /  dz  l  dx  dy  J 

(4.54) 


Once  more  exploiting  the  z-independence  of  the  susceptibility  components, 


vx(zM^d)  =  (“(gXZ22)-a(eXz‘2)) 


(4.55) 


)  _  (  Jlcpx  dX22  +  iky^gikox  _  gikpx  ^Xl2  \  , 

<9x  “  <9y  / 


Vx  j*  =  ft 


(4.56) 


Inserting  (4.48),  (4.51),  and  (4.56)  into  our  expression  for  the  first-order  response  electric 

E  E 

field  (4.46),  taking  Bo  —  —  — >  UqHq  =  —  — >■  Ho  =  ,  /—  Eq, 

c  c  V  Mo 


ErA  =  -<5  (to  -  cud) 


nEoi 

~Y~ 


.ft) 

%33  ~  l  — 
ft 


+  ik 0X22 


(4.57) 


But  k 


ft) 

— ,  so 
c 


-8  (cu-  coo) 


nEoi 

~Y~ 


I33  ~  ik 


dl22 

dx 


+  kok%22  +  ik^^J  zdrr 


(4.58) 


4.6.2  Second  Order  Response 

Returning  to  equation  (4.33),  the  second  order  response  electric  field  is 


£r,2=  I  Gk{r,r')E{  [  I  G,  (/ ,r")  F?  (e0,Ho)  d2r\ 


Gk(r',r")F2  (Hq,Eq)  d2r"  |  d 


i2-" 
(4.59) 
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The  first  integral  over  r  is  identical  in  structure  to  the  first  order  response  field,  given  in 
(4.58). 

J^Gk(r',r")F{"  (44)  d2r"  = 

-S(G)-c Oo)— ^  J  H()l)  (k\r'  -r"\)  elkox"  (k2X33  ~  +  kokX22  +  ik~§^J  zd2r" 

(4.60) 

Focusing  on  the  second  integral  over  r ' , 

J^Gk(r',r")F{"  (4,4)  d2r"  = 

Gk  (/,  r")  (-fc2  (x  *  4)  +  icoe o  (v  x  (xe  *  4)  )  -  V  (v  •  (xm  *  4)  )  )  d2/ 

(4.61) 

=  2nd  (co  -  (Oo)  J”  Gk  (?',  r")  eik^'  (-k2x*  ( -H0y )  +  lco£0V  x  (Xe  *  ( E0z )) 

-v(v-(r*(-^oy))))^" 

=  ^^<5  (CD  —  coo)  [  H^l)  (k\r-r"\)eikoX"  (k2,[^x*y  +  i(0£oyx{xe*z) 

L  ./—oo  y  y  jUo 

+i/?v(v-(^^»Vv' 

V  Mo  / 
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JZEo'i 


8  (ft)  -Wo)  f  H{Ql)  (k\r  -r'\)elk°x"  (k2<r^(Xi2X  +  X22y)  +  i(0£oV  x  (^33  z) 

J—oo  \  V  Mo 


+  \  ^^-(Xnx+Xny))  )  dzr 
V  Mo 


7r£’o/ 


8  (to 


coo)  (fclF'-r"!)^" 

J—oo 


—  {Xi2X  +  X22y)  +  i(0£o 


dX33  ,  8X33 


dy" 


dx 


-y 


/fbv  f  8xn  8x22 
VMo  V  dx  +  dy" 


o  // 

d2r 


Gk  (F/,r//)  (Ho,Eq)  d2r"  = 


%E()i 


5  (ft) 


-ft)0)  [  H^(k\r  ^r'\)elkoX"  (k2*r^(Xi2X  +  X22y)  +  iCD£o( 
j—oo  V  V  Mo  V 


<?Z33  ,  <?Z33 


5/ 


-  x  — 


dx" 


£o 

+  w— 

Mo 


^2Zl2  d~X22 


d(x")2  dx  dy 


x  + 


d2X  12  +  ^2Z22 


[d/dx"  d(y-y-_ 


y  dr  (4.62) 


Inserting  (4.60)  and  (4.62)  into  (4.59)  yields 
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Gk  (r,  r)  x 


n 


nE0i 


8((0- COq)7^-  I  //q1}  (k\r  -r"\)  e'k°x'  [k2X33  -  ik~^  +  k0kX22  +  zd2r" , 


5  (w 


Ob)  f  (*  I"  -  ?"|)  ^  (k2M  (Z12 

J —oo  y  y  jUQ 


—  (Zi2-?  +  X22;y)  +  z'a>£o 


<?Z33  ,  <?Z33  .  .  . 

X - + 


dy 


dx 


d2X\l  d2X22 


lE)\[d(x")  dx  dy 


x  + 


d2X 

dydx"  '  dyf 


12  <?"Z22 


d^f  |  d2r  (4.63) 


Er  2  = 


Gk  (r,  r)  8  (ft)  -  GJo)  \k2  KEq1 


x 


Z* 


1 1  (fc  |  /  -  r"  | )  e,fc()V"  (k2X33~ik  +  k0kX22  +  ik~^y  ^  zd2r"  )  )  + 


KEq 


(O/IqV  X 


m 


1}  (k\r-r"\)e  (Zi2*  +  Z22y)^  + 


i(0£o 

nE0i 


( dX33  <^Z33 


+ 


82Xn  d2X22 


dydx"  d(y"): 


y  + 


v(V  (x*  {k\r'  —  r"  | )  eik°x"  ^k2X33  ~  kk 


8X22 


dx  +k()kX22  +  ik-^^j 


/v  t2 — d  / 

zrf  r 


(4.64) 


We  may  observe  from  the  relative  complexity  of  the  above  result  that  the  second  and  higher 
order  terms  of  the  response  fields  (4.33)  and  (4.35)  will  be  computationally  taxing  for  all 
but  the  simplest  structures  and  material  parameters.  We  will  therefore  focus  our  efforts  in 
the  following  chapter  on  characterizing  first  order  responses  for  the  TO-derived  structures 
which  appear  in  the  literature. 
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CHAPTER  5: 
Applications 


In  this  chapter  we  will  apply  the  iterative  approach  outlined  in  Chapter  4  to  some  basic 
electromagnetic  redirection  structures  and  analyze  the  results  in  an  attempt  to  quantify  the 
error  fields  which  result  from  nonlinear  effects  in  these  TO  derived  devices. 


5.1  Cylindrical  Shield 

We  will  attempt  to  apply  the  iterative  method  of  Chapter  4  to  the  cylindrical  structure 
presented  in  [9]  in  Cartesian  coordinates,  followed  by  an  analysis  of  the  same  structure  in 
the  coordinate  system  most  suited  to  such  a  structure,  cylindrical  coordinates. 


5.1.1  Cartesian  Coordinates 

From  equation  (3.178)  the  electric  permittivity  £  is,  in  the  basis  (. x,y,z ), 


J 

/  .ft2  cos 

2  (0)  +  X2  sin2  (0) 

(■^2  —  pi)  cos  (0)  sin(0) 

0^ 

r 

£  =  £0 — 
Rr 

(ft2- 

pi)  cos  (0)  sin(0) 

ft2  sin2  (0)  +  yjcos2  (0) 

0 

V 

0 

0 

1/ 

(5.1) 


The  electric  susceptibility  Xe  corresponding  to  this  permittivity  is 


Xe 


1 

— £  - 
£o 


(l  0  0\ 
0  1  0 
Vo  0  l) 


(5.2) 


Xe  — 


^ ~r  COS2  (0)  +  ^7  sin2  (0)  -  1 

(^-^)cos(0)sin(0) 

V  0 


cos  (0) sin  (0) 
sin2  (0)  +  w  cos2  (0)  -  1 
0 


0  N 
0 


(5.3) 
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(5.4) 


(r^a  cqs2  (0)  +  _Z_  sin2  (0)  _  1  COS  (0)  sin  (0)  0 

Ze=  (™-^)cos(0)sin(0)  ^sin2(4>)  +  ^cos2(4>)-l 

V  0  0  if-1 

where  we  have  used  the  transformation  relation  (3.102)  in  the  last  step.  The  permittivity  is, 
as  for  the  case  of  the  square  structure  of  [8],  of  the  form 

(Xn  Xu  0  \ 

Ze  =  X21  X22  0  (5.5) 

\  o  0  X33 

and  so  our  result  (4.58)  for  the  first  order  response  to  an  incident  plane  wave,  polarized  in 
the  z-direction  will  apply  to  this  cylindrical  structure  as  well. 


Er  i  =  8  (<u  —  (Uq) 


7t  EqI 

~Y~ 


X33  ~  ik 


dl22 

dx 


+  kok%22  +  ik 


dXn 

dy 


^  zdrr' 
(5.6) 


Constructing  the  permittivity  as  per  [9]  (see  MATLAB  code,  Appendix  B)  for  a  structure  of 
inner  radius  1  meter  and  outer  radius  of  1 .5  meter,  we  note  immediately  from  Figure  5.1  the 
sharp  peaks  at  the  inner  radius  (set  to  1  meter  in  this  simulation).  The  mathematical  form  of 
the  electric  susceptibility  from  equation  (5.16)  suggests  that  these  peaks  are  an  artifact  of 
the  rectangular  grid  enforced  upon  an  inherently  cylindrical  structure.  That  the  location  and 
height  of  these  peaks,  in  addition,  varies  with  mesh  size  (see  Figure  5.1),  provides  further 
evidence  of  a  granularity  problem.  This  granularity  issue  and  the  resultant  inablility 
of  this  model  to  accurately  determine  susceptibility  values  according  to  the  TO 
prescription  suggests  that  any  application  of  the  iterative  approach  outlined  in  Chapter  4 
which  uses  these  susceptibility  values  should  be  suspect.  As  seen  in  Figure  5.2,  the 
attained  first  order  response  field  does  not  exhibit  any  of  the  expected  behavior  of  a  TO- 
derived  redirection  structure.  We  attribute  this  incongruence  to  the  issue  of  representing  a 
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a  circular  structure  on  a  rectangular  grid. 


X.r.r  for  Cylindrical  Redirection  Structure 


y  (meters)  -2  ,2  *  (meters) 


Xxx  for  Cylindrical  Redirection  Structure 


Figure  5.1:  Left:  xx  component  of  the  electric  susceptibility  as  computed  in  cartesian  coordinates, 
for  a  redirection  structure  of  inner  radius  1  meter,  outer  radius  1.5  meter,  and  mesh  size  of  0.05 
meter.  Right:  The  same  simulation,  with  a  mesh  size  of  0.03  meter.  Note  the  sharp  peaks  which 
arise  due  to  the  granularity  of  the  simulation  regardless  of  mesh  size. 


1st  Order  Response  Field  Amplitude  Relative  to  |2?o| 


1st  Order  Combined  Field  Amplitude  Relative  to  |2?o| 


Figure  5.2:  Left:  First  order  response  field  for  incident  £-polarized  plane  wave  with  wavenumber 
1  x  10  2 m  1 ,  in  Cartesian  coordinates,  mesh  size  0.03  meter  Right:  Combined  response  field. 
Note  that  the  jagged  peaks  of  Figure  5.1,  an  artifact  of  the  Cartesian  grid  imposed  on  the 
structure,  have  carried  over  to  the  response  field 


5.1.2  Cylindrical  Coordinates 

In  an  effort  to  address  the  granularity  issue  which  has  likely  compromised  the  usefulness 
of  the  treatment  seen  in  the  previous  subsection,  we  turn  to  the  coordinate  system  naturally 
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suited  to  a  cylindrical  structure  —  cylindrical  polar  coordinates.  Constructing  a  cylindrical 
coordinate  grid  and  inputting  the  susceptibility  values  determined  from  (3.189)  and  [9]  (see, 
once  again,  MATLAB  code,  Appendix  B),  we  observe  that  the  sharp  peaks  characteristic 
of  the  susceptibility  in  a  Cartesian  grid  system  appear  to  have  been  eliminated  in  favor  of 
the  expected  smoothly  varying  values  from  (3.178),  as  we  may  see  in  Figure  5.3. 


X.rj-  tor  Cylindrical  Redirection  Structure 

60 

40 

5  20 
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0 

-20 
2 

y  (meters)  -2  ,2  x  (meters) 

Figure  5.3:  xx  componenet  of  the  electric  susceptibility  for  the  same  cylindrical  structure  as  seen 

2k 

in  Figure  5.1,  in  cylindrical  coordinates.  Mesh  size  is  0.03  meter  in  the  radial  direction,  - 

6  200 
radian  azimuthally. 

From  equation(4.47),  for  an  incident  plane  wave  polarized  along  the  z-direction, 


E,-.\  =  <5(®-COo)^  f  Hq1)  (k\r-r\) 

A  J  —  oo 

(k2E0  (z  *  (VM'  z)  )  -  i (OiloHo  (V  x  (z  *  elkoX'  y)  )  +  V  (v  •  (z  *  (VV  z)  )  )  )  d2r' 

(5.7) 

Performing  the  computation  (5.9)  with  this  cylindrical  mesh,  we  note  that  although  the 
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overall  combined  (incident  and  first  order  response)  electric  field  amplitude  shows  overall 
promising  behavior,  falling  rapidly  within  the  structure,  the  resultant  fields  can  be  made  to 
vary  over  several  orders  of  magnitude  depending  upon  the  choice  of  mesh  fineness.  See 
Figures  5. 4-5. 6  for  a  comparison  of  the  combined  electric  field  for  various  mesh  sizes. 


1st  Order  Combined  Field  Amplitude  Relative  to  |J?q| 


y  (meters)  -2  « 

x  (meters) 


Figure  5.4:  Left:  Surface  plot  of  the  combined  (incident  and  first  order  response)  electric  field 

2k 

amplitude  for  the  cylindrical  structure.  Meshsize  is  0.003  meter  in  the  radial  direction, 

azimuthally.  Wavenumber  is  0.01  m-1.  Right:  The  response  field  for  the  same  meshsize,  along 
the  x-axis. 


1st  Order  Combined  Field  Amplitude  Relative  to  |£’o|  1st  Order  Combined  Field  Along  X-axis  Relative  to  |2£o| 

40  - 1 - i - 1 - i - i - . - i — 


Figure  5.5:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  for  the  same  structure  as  in 
Figure  5.4,  with  a  meshsize  of  0.005  meter  in  the  radial  direction. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 


900  ^ 
800  . 


Figure  5.6:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  for  the  same  structure  as  in 
Figure  5.4,  with  a  meshsize  of  0.007  meter  in  the  radial  direction. 


In  addition,  our  results  appear  strongly  dependent  upon  our  choice  of  incident  wavenumber, 
as  shown  in  Figures  5.7-5.10. 


m 


1st  Order  Combined  Field  Amplitude  Relative  to  |27o| 
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Figure  5.7:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  cylindrical 
redirection  structure  simulation,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  0.001m-1. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 


y  (meters)  -2  ,2  x  (meters) 


Figure  5.8:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  cylindrical 
redirection  structure  simulation,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  0.01m-1. 


1st  Order  Combined  Field  Amplitude  Relative  to  |i?o| 


y  (meters)  _2  ,2  *  (meters) 


1st  Order  Combined  Field  Along  X-axis  Relative  to  |i£o| 
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Figure  5.9:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  cylindrical 
redirection  structure  simulation,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  0.1m-1. 


We  note  from  Figure  5.10  that  as  the  wavelength  falls  to  less  than  three  orders  of  magni¬ 
tude  relative  to  the  radial  mesh  size,  the  simulation  appears  to  break  down  altogether.  We 
attribute  this  failure  to  the  inability  of  the  (now  relatively  coarse)  mesh  to  fully  extract  the 
necessary  data  from  our  model  for  the  incident  electromagnetic  wave. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 
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1st  Order  Combined  Field  Along  X-axis  Relative  to  | | 
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Figure  5.10:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  cylindrical 
redirection  structure  simulation,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  1/77  1 . 


5.2  Reduced  Material  Parameter  Cylindrical  Shield 


In  an  effort  to  eliminate  the  influence  of  diverging  material  parameters  upon  our  results, 

we  turn  to  an  approach  discussed  in  [3]  and  [11].  Enforcing  =  1,  Cummer  et  al. 

eij  uij 

have  shown  that  although  the  non-reflectance  condition  —  =  —  is  no  longer  met,  the 

£o  Mo 

remaining  relevant  material  parameters  in  cylindrical  coordinates  reduce  to 


(5.8) 


(5.9) 


The  absence  of  electrical  permittivity  components  other  than  £zz  will  prove  a  useful  sim¬ 
plification  when  considering  the  structure’s  response  to  a  plane  wave  polarized  along  the  z 


direction. 
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We  note  that  these  reduced  material  parameters  involve  only  one  cylindrical  coordinate 
component  which  varies  with  position  through  the  structure,  and  that  none  of  the  com¬ 
ponents  of  the  prescribed  susceptibilities  diverge  anywhere  within  the  structure  (see  Fig¬ 
ure  5.11).  These  reduced  parameters  are  therefore  expected  to  be  advantageous  from  a 
manufacturing  standpoint  [1 1],  as  arbitrarily  large  electromagnetic  material  parameters  are 
not  currently  realizable.  In  simulations,  the  reflected  intensity  arising  from  now  nonzero 
impedance  mismatch  between  surroundings  and  structure  is  somewhat  minor  [11]. 


Xxx  for  Cylindrical  Redirection  Structure 


Xyy  for  Cylindrical  Redirection  Structure 


* 


Figure  5.11:  Xm,xx  and  Xm,yy  for  the  reduced  material  parameter  cylindrical  structure.  Unlike 
the  material  parameter  solutions  of  the  previous  section,  these  susceptibiity  components  do  not 
diverge  as  we  approach  the  inner  boundary. 


Since,  for  this  application,  Xe  /  Xm,  the  general  first  order  response  (4.44)  becomes,  fol¬ 
lowing  development  identical  in  form  to  that  of  equations  (4.45)  to  (4.58), 


-8  (co-coq) 


nE0i 

~Y~ 


Xe,33 


ik^~  +  hkXm.22  +  ik^^m'2 


dx 


dy 
(5.10) 


Implementation  of  this  expression  for  the  first  order  response  electric  field  in  MATLAB 
(Appendix  B)  and  analysis  of  the  resulting  combined  (incident  and  first  order)  electric  field 
demonstrates  that  the  previously  observed  large  variation  of  results  with  respect  to  meshsize 
has  been  eliminated.  See  Figures  5.12-5.14  for  a  depiction  of  this  simulation  stability  with 


73 


respect  to  meshsize. 


1st  Order  Combined  Field  Amplitude  Relative  to  |£7q|  1st  Order  Combined  Field  Along  X-axis  Relative  to  |2£o| 
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Figure  5.12:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.05  meter,  azimuthal 
meshsize  SL  and  wavenumber  0.1  m-1. 
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Figure  5.13:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  O.lrn-1. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 
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Figure  5.14:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.008  meter,  azimuthal 
meshsize  and  wavenumber  0.1m-1. 


We  still  observe,  however,  a  very  strong  dependence  of  our  results  on  the  choice  of  incident 
wavenumber.  See  Figures  5.15-5.19  for  a  review  of  simulation  dependence  upon  incident 
wavenumber. 
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Figure  5.15:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  0.01m-1. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 
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Figure  5.16:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  0.1m-1. 
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Figure  5.17:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  lm-1. 


As  previously,  the  simulation  appears  to  breakdown  as  the  incident  wavelength  falls  to 
within  three  orders  of  magnitude  of  the  radial  meshsize. 


76 


1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 
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Figure  5.18:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  10m_1. 


1st  Order  Combined  Field  Amplitude  Relative  to  |f?o| 


y  (meters)  -2  ,2  x  (meters) 


Figure  5.19:  Surface  plot  and  two  dimensional  plot  along  the  x-axis  of  the  amplitude  of  the 
combined  (through  first  order)  electric  field,  relative  to  incident  field  amplitude,  for  the  MATLAB 
simulation  of  the  cylindrical  reduced  material  parameter  redirection  structure,  inner  radius  1 
meter,  outer  radius  1.5  meter,  with  cylindrical  mesh,  radial  meshsize  of  0.02  meter,  azimuthal 
meshsize  and  wavenumber  100m-1. 


Based  upon  the  above  results,  we  hypothesize  that  the  strong  dependence  on  meshsize  ob¬ 
served  in  Section  5.1  is  a  result  arising  from  the  divergence  of  the  material  parameters 
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observed  in  those  structures.  As  our  mesh  becomes  finer,  it  incorporates  a  greater  portion 
of  the  region  near  the  inner  boundary,  where  susceptibility  values  diverge,  leading  to  rising 
response  field  values,  as  we  have  seen.  The  simulation  results  obtained  in  the  case  of  the 
reduced  parameter  structure,  which  do  not  involve  diverging  electric  or  magnetic  suscepti¬ 
bility,  appear  highly  stable  with  respect  to  meshsize  variation.  A  common  feature  to  both 
approaches,  however,  is  the  strong  dependence  of  our  results  upon  incident  wavenumber. 


5.3  Square  Shield 

We  now  consider  the  square  electromagnetic  redirection  structure  outlined  in  [8].  Note  that 
as  for  the  unreduced  parameter  cylindrical  structure,  some  components  of  the  susceptibility 
diverge  as  we  approach  the  inner  surface.  In  Cartesian  coordinates,  the  susceptibility  tensor 
X  may  be  expressed  as 


(Xu 

X\2 

0  \ 

x  = 

X21 

Ill 

0 

0 

%33  / 

(5.11) 


For  an  incident  plane  wave  polarized  in  the  z-direction,  the  first  order  response  is  given  by 
equation  (4.58)  as 


Er  \  =S((0-C0o) 


H0l)  (■ k\r-r'\ ) 


Jkox' 


k2X33 


+  k0kX  22  + 


/v  j2->/ 

zd  r 


(5.12) 


Results  of  simulations  (Appendix  B)  for  this  structure  are  similar  to  those  observed  for 
the  Cylindrical  structure  once  an  appropriate  mesh  shape  was  imposed  upon  it.  We  again 
observe  similar  overall  behavior  between  all  structures,  displaying  some  of  the  desired 
combined  field  behavior  (field  amplitude  falls  to  a  minimum  rapidly  upon  entering  the 
inner  area  of  the  structure),  with  amplitudes  strongly  dependent  upon  meshsize,  as  shown 
in  Figures  5.20-5.22. 
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1st  Order  Combined  Field  Amplitude  Relative  to  |£7q| 


Figure  5.20:  Surface  plot  and  plot  along  the  x-axis  for  of  the  combined  (to  first  order)  electric 
field  amplitude,  relative  to  incident  field  amplitude,  for  the  square  structure  simulation,  with 
meshsize  0.03  meter,  wavenumber  0.02  nT  ,  half  side  lengths  1  and  2  meters. 
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Figure  5.21:  Surface  plot  and  plot  along  the  x-axis  for  of  the  combined  (to  first  order)  electric 
field  amplitude,  relative  to  incident  field  amplitude,  for  the  square  structure  simulation,  with 
meshsize  0.05  meter,  wavenumber  0.02  m  1 ,  half  side  lengths  1  and  2  meters. 
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1st  Order  Combined  Field  Amplitude  Relative  to  | Eq | 
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Figure  5.22:  Surface  plot  and  plot  along  the  x-axis  for  of  the  combined  (to  first  order)  electric 
field  amplitude,  relative  to  incident  field  amplitude,  for  the  square  structure  simulation,  with 
meshsize  0.07  meter,  wavenumber  0.02  m  1 ,  half  side  lengths  1  and  2  meters. 


In  addition  to  these  findings,  we  observe  combined  field  amplitudes  dependent  on 
wavenumber  as  per  the  cylindrical  structure.  Also  noted  is  the  breakdown  of  the  simu¬ 
lation  similar  to  that  of  Figures  5.10,  5.18,  and  5.19  upon  the  incident  wavelength  lowering 
to  within  three  orders  of  magnitude  of  the  meshsize. 
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CHAPTER  6: 
Conclusions 


Here  we  will  consider  the  results  of  the  approach  to  electromagnetic  scattering  pursued 
through  Chapters  4  and  5,  and  outline  possibilities  for  future  work. 


6.1  Effectiveness  of  the  Iterative  Approach 

Although  well-grounded  in  the  fundamental  relations  of  classical  electromagnetism,  im¬ 
plementation  of  a  Neumann  series  approximation  for  the  response  electric  and  magnetic 
fields  can  be  problematic  when  applied  to  exotic  structures  such  as  TO-based  CDEW  elec¬ 
tromagnetic  redirection  structures.  The  divergence  in  the  material  susceptibility  required 
for  even  the  simple  ideal  structures  we  have  examined  here  can  lead  to  considerable  vari¬ 
ation  in  the  results  of  our  numerical  calculations  depending  on  mesh  size  and  incident 
wavenumber.  For  this  reason,  it  is  likely  that,  in  order  to  attain  convergence  of  numeri¬ 
cal  results,  a  meshsize  finer  than  possible  on  the  hardware  employed  for  this  study  should 
be  used,  or  analyses  should  be  restricted  to  reduced  material  parameter  structures  such  as 
those  studied  in  Section  5.2.  Regardless,  powerful  hardware  would  also  be  required  to  be¬ 
gin  analyzing  the  numerical  results  for  second  and  higher  order  responses.  These  higher 
order  results  would  be  of  value  in  addressing  the  remaining  issues  of  series  convergence 
and  the  extreme  incident  wavenumber  dependence  of  our  results. 

Due  to  practical  engineering  concerns,  a  real-world  CDEW  redirection  structure  will  fea¬ 
ture  neither  a  continuously  varying  susceptibility,  nor  a  susceptibility  with  any  components 
diverging  as  a  particular  bounding  surface  is  approached.  Studies  such  as  [3]  and  [13]  have 
examined  the  linear  response  of  layered  structures,  designed  to  approximate  the  cylindrical 
redirection  device  proposed  in  [9]  and  analyzed  in  Chapter  3  above. 


6.2  Future  Work 

Opportunities  for  follow-on  research  are  plentiful.  As  mentioned,  attaining  a  convergent 
first  order  response  and  analyzing  the  numerical  computations  for  second  and  higher  order 
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response  fields  may  be  possible  with  access  to  greater  computing  power  than  was  available 
for  this  work.  In  addition,  once  the  higher  order  responses  have  been  well  characterized, 
and  assuming  that  the  Neumann  series  is  shown  to  converge,  the  efficacy  of  the  approach 
may  be  tested  via  experiment  using  the  layered  structures  which  have  been  shown  to  closely 
mimic  the  theorized  linear  behavior  of  TO-based  devices.  Such  experimental  work  could  be 
carried  out  as  an  extension  to  the  metamaterials  research  currently  being  done  by  Professors 
James  Luscombe  and  Dragoslav  Grbovic  of  the  NPS  Physics  Department. 

Beyond  these  refinements,  there  is  considerable  theoretical  ground  which  may  be  analyzed 
through  the  approach  taken  here.  In  order  to  be  considered  a  successful  extension  of  cur¬ 
rent  theory,  this  iterative  approach  should  allow  us  to  predict  observed  phenomena.  Miller’s 
rule,  an  empirical  method  of  attaining  higher  order  electric  susceptibilities  in  regular  non¬ 
magnetic  structures  (i.e.,  crystals)  from  the  lower  electric  susceptibilities,  has  been  shown 
to  be  effective  in  predicting  the  nonlinear  optical  behavior  of  these  structures.  If  valid, 
the  approach  taken  here  to  nonlinear  optical  responses  should  reduce  to  Miller’s  rule  in 
the  appropriate  material  regime.  As  a  further  test  of  the  predictive  power  of  our  approach, 
we  should  be  able  to  show  the  appearance  of  well-known  nonlinear  optical  effects  such  as 
frequency  doubling  and  the  Kerr  effect  in  appropriate  media. 
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APPENDIX  A: 
Green  Functions 


Green  functions  relevant  to  the  iterative  approach  of  chapter  5  of  this  work  will  be  derived 
here. 


A.l  The  Green  function  for  the  two  dimensional  Helmholtz 
equation 

The  homogenous  Helmholtz  equation  is  given  as 


V2/  +  ^/  =  0  (A.l) 

where  /  is  the  function  of  interest  and  ko  is  a  known  parameter.  The  Green  function  for  this 
equation  is  therefore  given  as  the  solution  to 

V2G  (r,  r ')  +kgG  (r,  rf)  =  8  (r-  r')  (A.2) 

Performing  a  Fourier  transform  and  applying  the  sifting  property  of  the  Dirac  delta  func¬ 
tion, 


(V2G(r,r')+k5G  (?,?')) 


~ik'7d2r  = 


(A. 3) 


e~lh7d2r+  [  k20G(r ,7) 

J  — oo 


ik?d2r  =  e 


ik-7J 


From  Green’s  second  identity. 


(A.4) 


[D{fV2g-gV2f)dA 


g{Vf-n))dl 


(A. 5) 


83 


where  D  is  the  domain  of  the  functions  /  and  g,  dl)  is  the  bounding  set  to  that  domain,  and 
n  is  the  normal  to  the  bounding  set.  Thus, 


[  fV2gdA  =  [  gV2fdA  —  [  (gd-l-fdJ-)dl 
Id  Jd  JdD  V  on  dn 


(A. 6) 


where  —  =  V/  •  n  and  similar  for  — .  Taking  /  =  e~,k  r  and  g  =  G,  and  furthermore  taking 
dn  dn 

D  as  the  entire  x-y  plane,  and  assuming  that  G  falls  off  at  least  as  fast  as  -,  assuming  a 

r 

value  of  zero  in  the  large  r  limit. 


e-il?V2G  (r,r)dA 


(A.7) 


e-il7V2G(r,r)dA 


(A. 8) 


e~il7V2G(r,r')dA 


(A. 9) 


The  two-dimensional  Fourier  transform  of  the  Laplacian  of  the  Green  function  is  thus  sim¬ 
ply  the  Fourier  transform  of  the  Green  function  multiplied  by  —k2,  where  k  is  the  Fourier 
transform  parameter.  Inserting  this  result  into  (A.7),  with  g{^  I-™  G  (G  ?')  e  lk  1  d2r. 


(, ki  —  k2)  G  (k,  r)  =  e~il7>  (A.  10) 

(A-U) 

The  Green  function  may  now  be  attained  by  inverse  Fourier  transforming  the  above  expres¬ 
sion. 
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(A.  12) 
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(A. 13) 


In  Cartesian  coordinates  in  k-space  (this  integral  may  also  be  performed  in  cylindrical 
coordinates), 
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where  p  =  r  —  r'. 


G 


1 

(2tt)2 


JkxPx 


k2 


(kl~ky) 


dkxdky 


(A. 14) 


(A. 15) 


G 


(A. 16) 


The  denominator  of  the  kx  integral  has  simple  poles  given  by  kx  —  ±  V  k2  —  k2.  If  kp  is 
greater  than  ky,  these  poles  lie  on  the  real  pv-axis.  If  kp  is  less  than  ky,  the  poles  lie  on  the 
imaginary  kx- axis.  If  kp  =  ky,  there  is  a  single  pole  of  second  order  at  the  origin.  We  will 
evaluate  (A.  16)  case  by  case. 


A.1.1  Case  1:  |&o|  >  |  ky  j 

The  integrand  of  the  kx  integral  has  poles  on  the  real  kx  line  at  kx  =  ±  ^  jk q  —  k2.  We  will 
evaluate  it  via  a  contour  in  the  complex  kx  plane  comprised  of  a  large  semicircle  in  the 
upper  half  of  the  plane,  taking  J ky}  —  kj  — »  Jk^—kj  +  i£, 
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Taking  f{z)  = 


,tzpx 


z  +  \/t$-lc*  +  ie  ]  \z- 


k q  —  ky  —  ie 


,  we  note  that  the  poles  of  f(z) 


now  lie  at  y/kf}  —  kj.  +  ie  and  —  v/ k^  —  k2  —  ie,  and  that  therefore  only  the  simple  pole 


Lo 


£q  —  k2  +  ie  lies  within  our  contour.  From  the  residue  theorem. 
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2m Res  (f(z  =  \/k2  -k2  +  ie  )  )  =  [  f  (kx)  dkx  +  /_  / (z)  dz 

-R 


(A. 17) 


where  Cr  denotes  the  semicircle  in  the  upper  half  plane  of  radius  R. 
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Taking  z  =  on  Cr, 
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As  R  — y  oo,  the  second  integral  approaches  a  value  of  zero  by  Jordan’s  lemma,  since 
sin(0)  >  0  on  0  =  [0,  n\. 
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Inserting  this  result  into  (A. 22), 
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A.1.2  Case  2:  |&o|  <  \ky\ 

In  this  case,  /  (z)  has  simple  poles  at  kx  —  ±i  y 

large  semicircle  in  the  upper  half  plane,  only  i 
residue  theorem, 

2 KiRes  (f  (z  =  i yjkj-k^j  )  =  J  ^  f  (kx)  dkx  +  f  (z)  dz  (A.30) 

where  f(z)  —  7 - ,  e'f  ' - ,  ,  and  Cr  is  again  the  semicircle  of  radius  R  in  the 

upper  half  plane,  centered  at  the  origin.  Taking  z  =  Re‘°  on  Cr, 


ky  —ky.  Again  closing  our  contour  with  a 
Iky —k^  lies  within  the  contour.  From  the 
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Inserting  this  result  into  (A. 34), 


(A. 37) 


(A. 38) 


(A. 39) 
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giPxk-x 


(A.40) 


(A.41) 


A.1.3  Case  3:  |&o|  =  \ky\ 

In  this  case,  our  kx  integral  reduces  to 


■jikxpx 
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(A.42) 


Moving  the  (now  second-order)  singularity  off  of  the  real  axis  as  for  case  1,  we  consider 


SZPx 
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where  C  is  once  more  the  large  semicircle  in  the  upper  half  plane.  By  the  residue  theorem. 
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where  /  (z)  = 
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As  previously,  the  integral  over  Cr  approaches  zero  as  R  — »  by  Jordan’s  lemma. 
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The  pole  at  z  —  is  is  simple,  and  therefore 
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(A.49) 


pikxPx 

_ 

(kx  +  is)(kx-is)  X 


2  ni  lim 


Z — 


■JPxZ 


[  (z  +  is) 


(A. 50) 


pikxPx 

_ Ab 

(kx  +  is)  (kx  -  is)  x 


2  ni 


eiPxi£ 

2  is 


(A. 51) 


pikxPx 

_ Ah 

(kx  +  is)(kx-is)  x 


ne  Px£ 


s 


(A. 52) 


This  integral  therefore  does  not  converge  as  s  — »  0.  We  conclude  that  ky  =  ko  is  an  unphys¬ 
ical  value  for  this  system. 


A.  1.4  Returning  to  the  Green  function 

From  (A.  16),  taking  all  allowed  values  of  ky  from  the  previous  three  subsections, 


(A. 53) 
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G(r,r) 


1  f  nielkypyepx^k°  ky  1  f  7ielkyPye  Px\iky  k° 

- ^  /  - .  - - dky  H - ^  /  - ,  - - dky 

(2?r)  %|<|*ol  (2?r)  %|>|*o|  ^/kj-kl 

(A. 54) 


1 

(2k)2 


■ikyPy+iPx^kl-kj 


dky 


+ 


(A. 55) 


(2nY 


“N  Ke  Pxvky  k°  (cos (. kypy )  +  i sin (. kypy )) 


_  t-2 
Kv  /CQ 


( 1 k  \ 


+ 


|*o I  ni  ( cos  ( kypy  +  px  J kl  -kj)+  i sin  ( kypy  +  px  J k%  -  tf 


J-\k0\ 


k-l  ~  kj 


dkv 


+ 


ne  Px\Jky  kQ  (cos  (kypy)  +  i  sin  (kypy)  ) 


%l 


ky  -  kl 


dky  ]  (A. 56) 


We  note  now  that  the  first  and  third  integrals  above  may  be  combined  to  yield  a  single 
integral  which  is  symmetric  about  the  origin,  while  the  second  integral  is  already  symmetric 
about  the  origin.  Further  splitting  these  integrals  into  even  and  odd  components. 
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(A. 57) 


where  Ti  is  the  entire  real  ky  line  except  for  the  interval  [—  |&o| ,  |&o|]>  which  is  symmetric 
about  the  origin.  Since  the  second  and  fourth  integrals  above  are  odd  in  ky  and  are  integrated 
symmetrically  with  respect  to  the  origin,  they  integrate  to  zero.  We  are  left  with 


G(r,/)  = 


ne 


Px \jky  k 


COS  (kyPy) 


(271)'- 


ky—ky 


dky  + 


M  nic os  (kypy  +  px  ^Jkl-kf) 

'-\k0\ 


dkv 


kl  -  kj 


(A. 58) 


The  remaining  integrals  are  each  even  and  integrated  symmetrically  with  respect  to  the 
origin. 


G 


2 

(2k)2 


Ke 


Px^kj  A'j 


COS  ( kypy J 


k2-kl 


dky  + 


"1*0  I  Ki  cos  kypy+p . 


k2-k2 


kl~k2 


(A. 59) 

Without  loss  of  generality,  let  p  lie  along  the  ky  axis.  Then  px  —  0  and  py  —  p  =  \r  —  r'\. 
We  are  left  with 
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G(r,r)  = 


(2  nY 


KC os  ( kyp ) 


ko 


I) 


:  dky  “h 


nic  os  (Xyp) 


\ 


*V*-( 


:  dk-v 


_  (  b  V 
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(A.60) 


K 

Let  w  =  -2.  Note  that  for  — >  Ifcol,  w  — >  ±1. 

k0 


/-» -»/\  1  A  [°°  7t cos  {k0pu)  f±l  mcos(k0pu 

G{r,r  >  =  2^2  l  /i  “ — /  „  .  *0«“  +  /  — — 7,  .  k()du 


koV  u~  —  1 


/o  koVl  —  w2 


(A. 61) 


Since  the  second  integral  above  is  even  in  u. 


f-y  (->  ~*t\  1  /  [°°  ncos(k0pu)  r1  nicos(k0pu 

Gy^r)  =  27^\Jl  — — /  n  .  k0du+  I  — — u  „  k0du 


ko  \/  u2  —  1 


'0  A:oV  1  —  M2 


(A. 62) 


For  the  second  integral  above,  let  u  —  cos  (/). 


G{?J)  =  — 


2tt  \  J l  Vw2  —  1 


cos(kopu)  ,  r° —i  cos  (kop  cost  t))  .  .  .  ,  ,  . 

v  K  J  du+  I  - ,  K  sin  (?)  dt  |  (A. 63) 


\/ sin2  (0 


„  _,x  1  (  [°°  cos  (kopu)  f  2 

G(r,r)  —  —  [  /  —  7  du  +  i  cos  (kop  cos  (t))  dt 

2 7T  V./l  Vm2  -  1  2o 


(A. 64) 


The  first  integral  above  is  an  integral  representation  of  the  Bessel  function  of  the  second 
kind,  order  0, 


f  COS  (k.Qpil)  7t  .  .  /zc\ 

/  /-  ^  du  =  -N0(k0p )  (A. 65) 

Vu2  - 1  2 

While  the  second  integral  of  (A. 63)  is  an  integral  representation  of  the  Bessel  function  of 
the  first  kind,  order  0, 
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(A. 66) 


r  2  ji 

/  cos  (kop  cos  (t))  dt  =  —Jo(kp) 

Jo  2 


Inserting  (A. 64)  and  (A. 65)  into  (A. 63), 


G(r,r')  =  ^A/o(M)  +  y-A)(*oP)^ 

G  (^)  =  \  Vo  (fc0 P)  Wo  (*Op)) 

G(aF0  =  ^dVop) 


where  H, 


(i) 


o 


is  the  first  Hankel  function  of  order  zero. 


(A. 67) 

(A. 68) 

(A. 69) 
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APPENDIX  B: 
MATLAB  Simulations 


This  appendix  will  present  the  MATLAB  scripts  used  to  simulate  the  various  TO-derived 
structures  studied  in  Chapter  5  of  this  work. 


B.l  Cylindrical  Cloak  in  Cartesian  Mesh 

%%  An  Iterative  Approach:  Cylindrical  Structure  with  Cartesian  Mesh 

o,  o. 
o  o 

%  Input  physical  constants 
%clearvars ; 

run ( ' Phy s ica ICons t ants . m' ) ; 

o,  o_ 
o  o 

o,  o. 
o  o 

%  Input  known  material  values 
a  =  1; 
b  =  1.5; 

R  =  (b-a) /b; 

o,  o_ 
o  o 

%  Input  meshgrid  for  the  region  of  interest, 
increment  =  0.03; 

x  =  cat (2, -b : increment : -increment,  increment : increment :b)  ; 
y  =  cat (2, -b : increment : -increment, increment : increment :b) ; 

[X,Y]  =  meshgrid (x,  y)  ; 

o.  o_ 
o  o 

%  Convert  to  cylindrical  coordinates 

r  =  sqrt  (X.A2+Y.A2)  ; 

origin  =  size(x,2)/2; 

theta  =  zeros ( size (x, 2 ) , size (x, 2 ) ) ; 

for  n  =  1 : 1 : size  (x, 2 ) /2 

for  m  =  1 : 1 : size (x, 2) /2 

theta (n,m)  =  atan(y(n)/x(m))+pi; 

theta (n+size (x, 2 ) /2 , m)  =  atan (y (n+size (x, 2 ) /2 ) /x (m) ) +pi; 
theta (n, m+size (x, 2 ) /2 )  =  atan (y (n) /x (m+size (x, 2 ) /2 ) ) +2 *pi ; 
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theta (n+size (x, 2 ) /2 , m+size (x,  2 ) /2 )  =  ... 

atan (y (n+size (x, 2 ) /  2) /x (m+size (x,  2 )  / 2)  )  ; 


end 


end 


%  Input  the  coordinate  transformation  from  $r$  to  $r/' 
r_prime  =  (r-a)  . /R; 
for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size (x, 2 ) 
if  r (n, m)  <=  a 

r_prime(n,m)  =  0; 

end 

end 

end 


{  f  }$■ 


%  Input  the  transformation  optics  derived  susceptibilities.  Assume  vacuum 
%  outside  of  the  "shield." 

L_1 1  =  zeros (size (x, 2 ) ,  size (x,  2 ) )  ; 

L_12  =  zeros (size (x, 2 ) ,  size (x, 2 ) ) ; 

L_13  =  zeros (size (x,  2 )  ,  size (x,  2 )  )  ; 

L_22  =  zeros ( size (x, 2 ) ,  size  (x, 2 ) ) ; 

L_23  =  zeros (size (x,  2 )  ,  size (x,  2 )  )  ; 

L_33  =  ones (size (x, 2 ) ,  size (x, 2 ) ) ; 

L  =  cell (size (x, 2 ) ,  size (x,  2 ) )  ; 

Chi_ll  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_12  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_13  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_22  =  zeros (size (x, 2) , size (x, 2) )  ; 

Chi_21  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_23  =  zeros ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_31  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_32  =  zeros  ( size (x, 2 ) , size (x, 2 ) ) ; 

Chi_33  =  zeros (size (x,2) , size (x, 2) ) ; 
for  n  =  l:l:size(x,2); 

for  m  =  1 : 1 : size  (x, 2 ) 

L_ll(n,m)  =  R+ (a/r_prime (n^ m) )*( sin (theta (n, m) )) A2 ; 

L_12(n,m)  =  -a*cos (theta (n, m) ) *sin (theta (n, m) ) /r_prime (n, m) ; 
L_22(n,m)  =  R+ (a/r_prime (n^ m) )* (cos (theta (n, m) )) A2 ; 

end 
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86 
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end 

for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size  (x, 2 ) 

L{n,m}  =  [L_ll (n,m) , L_12 (n, m) , L_13 (n,m) ; L_12 (n, m) , L_22 (n,m) , . . . 
L_23  (n,  m)  ;  L_13  (n,  m)  ,  L_23  (n,  m)  ,  L_33  (n,  m)  ]  ; 

end 


end 


Chi  =  cell (size (x, 2 ) ,  size (x, 2 ) ) ; 
for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size  (x, 2 ) 
if  r (n, m)  >  b 

Chi{n,m}  =  zeros (3,3); 
else  if  r(n,m)  <=  a 

Chi{n,m}  =  zeros(3,3); 

else 

Chi{n,m}  =  (L { n, m} * ( (L { n, m} ) 
r  (n,  m)  )  )  -eye  (3)  ; 

end 

end 

end 


end 

f or  n  = 
for 


1:1: size  (x, 2 ) 
m  =  1 : 1  :  size  (. 
Chi_l 1 (n, m)  = 
Chi_12 (n, m)  = 
Chi_13 (n, m)  = 
Chi_2 1 (n, m)  = 
Chi_22 (n, m)  = 
Chi_23 (n, m)  = 
Chi_31 (n, m)  = 
Chi_32 (n, m)  = 
Chi_33 (n, m)  = 


,2) 

Chi { n, m} (1,1) 
Chi { n, m} (1,2) 
Chi { n, m} (1,3) 
Chi { n, m} (2,1) 
Chi { n, m} (2,2) 
Chi { n, m} (2,3) 
Chi { n, m} (3,1) 
Chi { n, m} (3,2) 
Chi { n, m} (3,3) 


end 


) ) . * ( r_prime (n,m) / (R*. . . 


end 

o,  o. 
o  o 

%  Input  the  incident  fields,  as  well  as  the  functions  $F_1 (E_0, H_0 ) $  and 
%  $F_2  (H_0 ,  E_0 )  $ .  Begin  by  choosing  a  value  for  $\omega_0$.  Take  $10/'{11}  Hz$. 
k_0  =  le-2; 


o,  o, 
o  o 
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%  Input  the  initial  (frequency  domain)  electric  field 
E0_z  =  2 *pi . *exp ( li*k_0 . *X) ; 

o.  o_ 
o  o 

%  Compute  the  tensor  contractions  appearing  in  $F_1$  and  $F_2$ 


o,  o. 
o  o 

%  Perform  the  integration  indicated  in  equation  (3.40)  to  obtain  the 
%  first-order  response  electric  field.  Begin  on  the  x-axis  outside  the 
%  shield.  That  is,  let  $\vec{r}  =  -2a\hat{x}$. 
rminusrprime  =  cell (size (x, 2) , size  (x, 2) ) ; 

G1  =  cell ( size (x, 2 ) , size (x, 2 ) ) ; 

51  =  3*X.A2.*Y.A4+3*X.A4.*Y.A2-r.A5+X.A6+Y.A6; 

52  =  X. A2-2*r . A3+Y. A2; 

Chi_22_x  =  X.* (6*X . A2 . *r-r . A3+2*Y. A2) ./Sl-X.* (12*X. A2 . *Y. A2-5*r . A3+6* . . . 

X . A4+6*Y . A4 )  . * (X. A2+Y. A2-r. A5+2 . *X . A2 . *r. A3+Y. A4)  ./  (S1.A2) ; 

Chi_12_y  =  X. *S2 . /S1-2*X. *Y. A2 . * (3*r-l) . /Sl-X . *Y . A2 . *S2 . * ( 12*X . A2 . *Y . A2- . . . 
5*r . A3+6*X . A4+6*Y. A4) ./(S1.A2); 


for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size (x, 2 ) 
if  r (n, m)  >  b 

Chi_22_x (n, m)  =  0; 

Chi_12_y (n,  m)  =  0; 
elseif  r(n,m)  <=  a 

Chi_22_x (n,  m)  =  0; 

Chi_12_y (n,  m)  =  0; 

end 

end 

end 

for  n  =  1 : 1 : size  (x, 2 ) ; 

for  m  =  1 : 1 : size (x, 2 ) ; 

rminusrprime { n, m}  =  sqrt ( (X (n, m) -X) ,A2+(Y(n,m)-Y) .  A2 )  ; 
Gl{n,m}  =  besselh (0, 1, k_0 . *rminusrprime{n,m} )  ; 

G1  { n,  m}  (n,  m)  =  0  ; 
end 

end 

o,  o. 
o  o 

%  Compute  the  integral  (4.58) 

F  =  cell(size(x,2),size(x,2)); 
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for  n  =  1 : 1 : size  (x, 2 ) 

for  m  =  1 : 1 : size (x, 2 ) 

F{n,m}  =  pi*li/2*Gl {n,m} . *exp (li*k_0*X) . * (k_0 A2*Chi_33-li*k_0*  ... 
Chi_22_x+k_0/'2 *Chi_22  +  li*k_0*Chi_12_y )  ; 

end 

end 

E0_z  =  2*pi*exp ( li*k_0*X)  ; 

E0_x  =  0; 

El_z  =  zeros (size (x, 2 ) ,  size (x, 2 ) ) ; 
for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size  (x, 2 ) 

El_z(n,m)  =  increment A2*trapz (trapz (F { n, m} ,  2 )) ; 

end 

end 

o,  o. 
o  o 

%  Linear  approach. 

E_0  =  1; 

D_x  =  (Chi_ll+1) *E_0; 

D_y  =  (Chi_21) *E_0; 

D_r  =  D_x  . *cos (theta) +D_y . *sin (theta)  ; 
for  n  =  1 : 1 : size (x, 2 ) 

for  m  =  1 : 1 : size (x, 2 ) 
if  r (n, m)  <=a 

D_r (n, m)  =  0 ; 

end 

end 

end 

o,  o, 
o  o 

%  Streamline  plot  the  linear  field 
streamline (X, Y, D_x, D_y, -b*ones (l,size(x,2)),y); 

o,  o. 
o  o 

%  Graphs 

hi  =  surf (X, Y, Chi_l 1 ) ; 

title (' $\chi_{xx} $  for  Cylindrical  Redirection  Structure interpreter . 
'latex ' , ' font size ' , 14 ) ; 

xlabel ( ' x  (meters )  ' ,  ' interpreter '  ,  'latex '  ,  ' font size '  ,  12 )  ; 
y label ( ' y  (meters )  '  ,  ' interpreter ' ,  'latex ' ,  ' font size ' , 12 ) ; 
z label ( ' $\chi_{ xx } $ ! , ' interpreter ' , 'latex ' , ' font size ' , 14 ) ; 
figure 
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surf (X, Y, abs (E0_z+El_z) / (2*pi)  ) 

xlabel ( ’ x  (meters )  ' ,  ' interpreter ' ,  ' latex '  )  ; 

ylabel ( ' y  (meters )  ' ,  ' interpreter  ' ,  'latex ') ; 

z label ( ' $\frac{\left | E_0+E_1 \ right | } { \left | E_0\ right | }$ 1 ,  ' interpreter ' , . . . 

'latex ' ,  ' font size ' , 14 ) ; 
figure 

surf (X, Y, abs (El_z) / (2*pi)  )  ; 

title ('1st  Order  Response  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' ,  . . . 

' interpreter ' ,  ' latex ' ) ; 
xlabel ( ' x  (meters ) ' ,  ' interpreter ' ,  ' latex ' ) ; 
ylabel ( ' y  (meters ) ' , ' interpreter ' ,  'latex ') ; 

z label ( ' $\f rac { \left | E_l\ right | } { \left | E_0\ right | } $ ' , ' interpreter ' ,  . . . 

'latex ' ,  ' font size ' , 14 ) ; 


B.2  Cylindrical  Structure  in  Cylindrical  Mesh 

%%  An  Iterative  Approach:  Cylindrical  Structure  in  Cylindrical  Mesh 

o,  o. 
o  o 

%  Input  physical  constants 
clearvars ; 

%run ( ' PhysicalConstants . m ' ) ; 

2-  2- 
o  o 

%  Input  known  material  parameters 
a  =  1; 
b  =  1.5; 

R  =  (b-a) /b; 

2-  2- 
o  o 

%  Input  cylindrical  meshgrid  for  the  region  of  interest 
rstep  =  0.02; 
rmax  =  2 ; 
phistep  =  pi/50; 

[r,phi]  =  meshgrid ( 0 : rstep : rmax, 0 : phistep : 2*pi)  ; 

o.  o. 
o  o 

%  Input  relative  permittivity  in  cylindrical  coordinates. 

eps_rr  =  ( (r-a) . /r) . A2 ; 

eps_pp  =  ones (size (r) ) ; 

eps_zz  =  (1/R) A2 . *ones (size (r) ) ; 
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n  =  1: 

1 :  size (phi , 1 ) 

for  m 

=  1:1: size  (r, 2) 

if 

r (n,  m) >b 

eps_rr (n, m)  = 

i; 

eps_pp(n,m)  = 

i; 

eps_zz (n, m)  = 

i; 

el 

seif  r  (n, m) <=a 

eps_rr (n, m)  = 

i; 

eps_pp(n,m)  = 

i; 

eps_zz (n, m)  = 

i; 

end 

end 


end 

o,  o. 
o  o 


%  Convert 

eps_xx  = 
eps_yy  = 
eps_xy  = 
chi_xx  = 
chi_yy  = 
chi_xy  = 
chi_yx  = 
chi_zz  = 

o,  o. 
o  o 


to  Cartesian  coordinates.  Compute  the 
eps_rr  .  *  ( cos  (phi )  )  .  A2+eps_pp  .  *  ( sin  (phi ) 
eps_rr . * ( sin (phi ) ) . A2+eps_pp . * (cos (phi ) 
(eps_rr-eps_pp)  . *sin (phi)  . *cos (phi)  ; 
eps_xx-l ; 
eps_yy-l ; 
eps_xy ; 
eps_xy; 
eps_z  z  — 1 ; 


susceptibility . 
)  .  A2; 

)  .  A2; 


%  Plot  the  xx  component  of  the  permittivity  for  comparison  with  the  purely 
%  Cartesian  approach 

hi  =  surf (r . *cos (phi) , r . *sin (phi) , chi_xx) ; 
hold  on 

title (' $\chi_{xx} $  for  Cylindrical  Redirection  Structure interpreter . 
'latex ' , ' font size ' , 14 ) ; 

xlabel ( ' x  (meters )  ' ,  ' interpreter '  ,  'latex '  ,  ' font size '  ,  12 )  ; 
y label ( ' y  (meters)',’ interpreter ' , 'latex ' , ' font size ' , 12 ) ; 
z label ( ' $ \chi_{ xx } $ ' , ' interpreter ' , 'latex ' , ' font size ' , 14 ) ; 
hold  off 


o,  o, 
o  o 


%  Input  incident  wavenumber  and  electric  field 
k_0  =  le-1; 

E0_z  =  2 *pi*exp ( li*k_0 *r . *cos (phi ) )  ; 

o,  o. 
o  o 
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%  Input  the  terms  to  be  integrated 
X  =  r . *cos (phi ) ; 

Y  =  r . *sin (phi ) ; 

51  =  3*X. A2 . *Y. A4+3*X. A4 . *Y. A2-r . A5+X. A6+Y. A6; 

52  =  X.  A2-2*r . A3+Y. A2; 

chi_yy_x  =  X. * (6*X. A2 . *r-r . A3+2 *Y . A2 )  . /Sl-X. * (12*X. A2 . *Y . A2-5*r . A3  +  6* .  .  . 

X. A4  +  6*Y.A4)  . * (X. A2+Y. A2-r . A5+2 . *X. A2 . *r . A3+Y. A4)  . / (SI . A2)  ; 
chi_xy_y  =  X. *S2 . /S1-2*X . *Y . A2 . * (3*r-l)  . /Sl-X. *Y. A2 . *S2 . * ( 12*X . A2 . * . . . 

Y. A2-5*r.A3+6*X.A4+6*Y.A4)  . /  (SI . A2) ; 
for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size  (r, 2 ) 
if  r (n, m)  >  b 

chi_yy_x (n, m)  =  0; 
chi_xy_y (n, m)  =  0; 
elseif  r(n,m)  <=  a 

chi_yy_x (n, m)  =  0; 
chi_xy_y (n, m)  =  0; 

end 

end 

end 

rminusrprime  =  cell (size (r) ) ; 

G  =  cell (size (r) ) ; 
rdphi  =  phistep*r ( 1 ,  : )  ; 
for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size  (r, 2 ) 

rminusrprime { n, m}  =  sqrt ( (X (n, m) -X) . A2+ ( Y (n, m) -Y) . A2 ) ; 

G{n,m}  =  -besselh (0, 1 , k_0 . *rminusrprime { n, m} ) ; 

G{n, m} (n, m)  =  0; 

end 

end 

F  =  cell (size  (r) ) ; 
for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size  (r, 2 ) 

F{n,m}  =  pi*li/2*G{n, m} . *exp (li*k_0*X) . * (k_0 A2*chi_zz-li*k_0*  .. 
chi_yy_x+k_0  A 2  * chi_yy  + 1  i * k_0  * chi_xy_y )  ; 

end 

end 

El_z  =  zeros (size (r) ) ; 
for  n  =  1 : 1 : size  (r , 1 ) 
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for  m  =  1 : 1 : size  (r, 2 ) 

El_z(n,m)  =  rstep*trapz ( rdphi . *trapz (F { n, m} , 1 ) ) ; 

end 

end 

o.  o. 
o  o 

%  Graphs 
figure; 

h2  =  surf (X, Y, abs (El_z+E0_z) / (2*pi) ) ; 
hold  on 

title ('1st  Order  Response  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' , . . . 

' interpreter latex ') ; 
xlabel ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 
y label ( ' y  (meters ) ' , ' interpreter ' , 'latex ' ) ; 

z label ( ' $\frac{\left | E_0+E_1\ right | } { \left | E_0\ right | } $ ' , 'interpreter ' , . . . 
'latex ' , ' font size ' , 14 ) ; 

title ('1st  Order  Combined  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' , . . . 

' interpreter ' , ' latex ' ) ; 
hold  off 

xx  =  -rmax : rstep : rmax; 

E0_c  =  cat (2, fliplr (E0_z (51,  : ) )  ,  E0_z (1,  : ) ) ; 

El_c  =  cat (2, fliplr (El_z (51, : ) ) , El_z (1, : ) ) ; 

E0_c (size (r,  2) +1)  =  []; 

El_c (size (r,  2) +1)  =  []; 
figure 

plot (xx, abs (E0_c+El_c) / (2*pi) ) 
hold  on 

T  =  '1st  Order  Combined  Field  Along  X-axis  Relative  to  $\left | E_0\right | $ ' ; 

title (T,  ' Interpreter ' ,  'latex ' ) ; 

xlabel ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 

ylabel ( ' $\frac{ \left | E_0+E_l\right |  } { \left | E_0\right I  } $ ' ,  ' interpreter '  ,  .  .  . 

'latex ' , ' font size ' , 14 ) ; 
hold  off 


B.3  Reduced  Material  Parameter  Cylindrical  Structure 
in  Cylindrical  Mesh 

%%  An  Iterative  Approach:  Reduced  Parameter  Cylindrical  Structure 
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o,  o. 
o  o 

%  Input  physical  constants 
clearvars; 

%run ( ' Phys i ca ICons t ant s . m 1 ) ; 

o.  o. 
o  o 

%  Input  known  material  parameters 
a  =  1; 
b  =  1.5; 

R  =  (b-a) /b; 

o,  o. 
o  o 

%  Input  cylindrical  meshgrid  for  the  region  of  interest 
rstep  =  0.008; 
rmax  =  2 ; 
phistep  =  pi/50; 

[r,phi]  =  meshgrid ( 0 : rstep : rmax, 0 : phistep : 2*pi)  ; 

o.  o. 
o  o 

%  Input  relative  permittivity  in  cylindrical  coordinates. 

eps_rr  =  ( (r-a) . /r) . A2; 

eps_pp  =  ones (size (r) ) ; 

eps_zz  =  (1/R) A2 . *ones (size (r) ) ; 


siz 

e (phi. 

1) 

1 : 1 

: size  ( 

r,  2) 

'  (n. 

m)  >b 

eps. 

_rr (n. 

m)  = 

i; 

eps. 

_PP  (n. 

m)  = 

i; 

eps. 

_zz (n. 

m)  = 

i; 

:  i  f 

r  (n,  m) 

<=a 

eps. 

_rr (n. 

m)  = 

i; 

eps. 

_PP  (n. 

m)  = 

i; 

eps. 

_zz (n. 

m)  = 

i; 

end 

end 

end 

o,  o. 
o  o 

%  Convert  to  Cartesian  coordinates.  Compute  the  susceptibility. 
eps_xx  =  eps_rr  .  *  ( cos  (phi )  )  .  A2+eps_pp  .  *  ( sin  (phi )  )  .  ^2  ; 
eps_yy  =  eps_rr . * ( sin (phi ) ) . A2+eps_pp . * ( cos (phi ) ) . ^2 ; 
eps_xy  =  (eps_rr-eps_pp)  . *sin (phi)  . *cos (phi)  ; 
chi_xx  =  eps_xx-l; 
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chi_yy  =  eps_yy-l; 
chi_xy  =  eps_xy; 
chi_yx  =  eps_xy; 
chi_zz  =  eps_zz-l; 

o,  o. 
o  o 

%  Plot  the  xx  component  of  the  permittivity  for  comparison  with  the  purely 
%  Cartesian  approach 

hi  =  surf (r . *cos (phi) , r . *sin (phi) , chi_xx) ; 
hold  on 

title (' $\chi_{xx} $  for  Cylindrical  Redirection  Structure interpreter . 
'latex '  ,  ' font size ' , 14 ) ; 

xlabel ( ' x  (meters )  '  ,  ' interpreter ' ,  'latex ' ,  ' font size ' , 12 ) ; 
y label ( ' y  (meters)',' interpreter '  ,  'latex ' ,  ' font size ' , 12 ) ; 
z label ( ' $\chi_{ xx } $  1 ,  ' interpreter '  ,  'latex '  ,  ' font size '  ,  14 )  ; 
hold  off 

o,  o. 
o  o 

%  Input  incident  wavenumber  and  electric  field 
k_0  =  le-1  ; 

E0_z  =  2 *pi*exp ( li*k_0 *r . *cos (phi ) )  ; 

o,  o. 
o  o 

%  Input  the  terms  to  be  integrated 
X  =  r . *cos (phi ) ; 

Y  =  r . *sin (phi) ; 

chi_yy_x  =  (2*X . *Y . A2 )  .*  (3*r-2)  ./  (r .  A 6) ; 

chi_xy_y  =  -X . * ( 2*X . A2 . *r-4 *Y . A2 . *r-X . A2+3*Y . A2 ) . / ( r . A 6 ) ; 
for  n  =  1 : 1 : size (r, 1 ) 

for  m  =  1 : 1 : size  (r, 2) 
if  r (n, m)  >  b 

chi_yy_x (n, m)  =  0; 
chi_xy_y (n, m)  =  0; 
elseif  r(n,m)  <=  a 

chi_yy_x (n, m)  =  0; 
chi_xy_y (n, m)  =  0; 

end 

end 

end 

rminusrprime  =  cell (size (r)  )  ; 

G  =  cell (size  (r)); 
rdphi  =  phistep*r  ( 1 ,  :  ) ; 
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for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size  (r, 2 ) 

rminusrprime { n, m}  =  sqrt ( (X (n, m) -X) . A2+ ( Y (n, m) -Y) . ^2 ) ; 

G{n,m}  =  -besselh ( 0 , 1 , k_0 . *rminusrprime { n, m} ) ; 

G{n,m} (n,m)  =  0; 

end 

end 

F  -  cell (size  (r) ) ; 
for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size  (r, 2 ) 

F{n,m}  =  pi*li/2*G{n, m} . *exp (li*k_0*X) . * (k_0 A2*chi_zz-li*k_0*  ... 
chi_yy_x+k_0 A2 *chi_yy+li*k_0*chi_xy_y )  ; 

end 

end 

El_z  =  zeros ( size  ( r) ) ; 
for  n  =  1 : 1 : size  (r , 1 ) 

for  m  =  1 : 1 : size (r, 2 ) 

El_z(n,m)  =  rstep*trapz (rdphi . *trapz (F { n, m} , 1 ) ) ; 

end 

end 

o,  o, 
o  o 

%  Graphs 
figure; 

h2  =  surf (X, Y, abs (El_z+E0_z) / (2*pi) ) ; 
hold  on 

title ('1st  Order  Combined  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' , . . . 

' interpreter ' , ' latex ' ) ; 
xlabel ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 
y label ( ' y  (meters ) ' , ' interpreter ' , 'latex ' ) ; 

zlabel ( ' $\frac{ \left | E_0+E_l\right | } { \left | E_0\right I } $ ' , ' interpreter ' , . . . 
'latex '  ,  ' font size ' , 14 ) ; 

title ('1st  Order  Combined  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' , . . . 

'  interpreter ' ,  ' latex ' ) ; 
hold  off 

xx  =  -rmax : rstep : rmax; 

E0_c  =  cat (2, fliplr (E0_z (51, : ) ) , E0_z (1, : ) ) ; 

El_c  =  cat (2, fliplr (El_z (51, : ) ) ,El_z (1, : ) ) ; 

E0_c (size (r,  2) +1)  =  []; 

El_c (size (r,  2) +1)  =  []; 
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figure 

plot (xx, abs (E0_c+El_c) /  (2*pi) ) 
hold  on 

T  =  '1st  Order  Combined  Field  Along  X-axis  Relative  to  $\left | E_0\right | $ ’ ; 

title  (T,  ' Interpreter ' ,  ’  latex ' ) ; 

xlabel ( ' x  (meters )  ’ ,  1  interpreter  1 ,  1  latex ' ) ; 

ylabel ( ' $\frac{\left | E_0+E_1\ right | } { \left | E_0\ right | } $ ' ,  ’ interpreter 1 ,  . . . 

’latex ' ,  ' font size ' , 14) ; 
hold  off 


B.4  Square  Cloak 


%%  An  Iterative  Approach:  Symbolic  Method  for  Square  Cloak 

o,  o. 
o  o 

%  Input  physical  constants 
clearvars ; 

run ( ' Physica ICon st ants . m ' ) ; 

g,  g. 
o  o 

%  Input  predetermined  structure  shape.  Let  $S_1$  and  $S_2$  be  the  inner  and 
%  outer  half  side  lengths,  respectively, 
si  =  1; 
s  2  =  2 ; 

g,  g, 
o  o 

%  Input  symbolic  expression  for  the  electric  permittivity,  as  per  Pendry  et 
%  al .  Rotate  the  symbolic  expression  through  $\f rac { \pi } { 2 } $ ,  $\pi$,  and 
%  $\frac{ 3\pi } { 2 } $  in  order  to  obtain  the  permittivity  throughout  the 
%  structure, 
syms  x  y  z; 
a  =  s2/ (s2-sl ) ; 

b  =  symfun  (a*y*sl/  (x/'2)  ,  [x  y]  )  ; 
c  =  symfun (a* (1-sl/x) , [x  y] ) ; 

epsl  =  symfun ([c/a,-b/a,0;-b/a, (aA2+bA2)/ (a*c),0;0,0,a*c], [x  y]); 

2-  9- 
o  o 

%  Input  meshgrid  for  the  region  of  interest, 
increment  =  0.03; 

[X, Y]  =  meshgrid (-1 . 5*s2 : increment :1.5*s2,-1.5*s2: increment : 1 . 5*s2 ) ; 

2-  2- 
o  o 
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%  Determine  numeric  values  of  the  permittivity  over  one  quarter  of  the 
%  structure,  and  rotate,  as  per  [8] . 

A_90  =  [0, -1,0; 1,0,0; 0,0,1] ; 

A_n90  =  [0,1,0; -1,0,0;  0,0,1]; 

A_180  =  [-1,0,0;0,-1,0;0,0,1]  ; 

A_n 18  0  =  A_1 8  0 ; 

A_2 7 0  =  [0, 1,  0;-l, 0, 0; 0,  0,  1]  ; 

A_n270  =  [0,-1, 0; 1, 0, 0; 0, 0, 1] ; 
eps2  =  symfun(A_90*  . . . 

eps 1 ( dot ( A_n  9  0  * [ x ; y ; 0 ] ,  [ 1 ; 0 ; 0 ] ) , dot ( A_n  9  0  * [ x ; y ; 0 ] ,  [ 0 ; 1 ; 0 ] ) ) * . . . 

A_90 . ' , [x  y]); 
eps3  =  symfun (A_180*  ... 

epsl (dot (A_nl80* [x; y; 0] , [1; 0; 0] ) , dot (A_nl80* [x;y; 0] , [0; 1; 0] ) ) *A_180 .',... 
[x  y]  )  ; 

eps4  =  symfun (A_270*  . . . 

eps 1 ( dot ( A_n2 7  0  *  [  x ;  y ;  0  ]  ,  [ 1 ; 0 ; 0 ] ) , dot ( A_n2 7 0 * [ x ; y ; 0 ] ,  [ 0 ; 1 ; 0 ] ) ) *A_2 7  0.',... 
[x  y]  )  ; 

chil  =  symfun (epsl (x, y) -eye (3)  ,  [x  y] )  ; 
chi2  =  symfun (eps2 (x, y) -eye ( 3 ), [x  y] ) ; 
chi3  =  symfun (eps3 (x, y) -eye ( 3 ),  [x  y] )  ; 
chi4  =  symfun (eps4 (x, y) -eye ( 3 ),  [x  y] )  ; 

CHI  =  cell (size (X, 2) , size (X, 2) ) ; 

CHI_1 1  =  zeros  ( size (X, 2 ) , size (X, 2 ) ) ; 

CHI_12  =  zeros  ( size (X, 2 ) , size (X, 2 ) ) ; 

CHI_22  =  zeros  ( size (X, 2 ) , size (X, 2 ) ) ; 

CHI_33  =  zeros  ( size (X, 2 ) , size (X, 2 ) ) ; 

G  -  cell (size (X, 2) , size  (X, 2) ) ; 

k_0  =  2e-l ; 

for  n  =  1 : 1 : size  (X, 2 ) 

for  m  =  1 : 1 : size (X, 2 ) 

if  atan2 ( Y (n, m) , X (n, m) )  >=  -pi/4  &&  atan2 ( Y (n, m) , X (n, m) )  <  pi/4 
if  abs(X(n,m))  ==  si 

CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI{n,m}  =  zeros(3,3); 

else 

CHI{n,m}  =  double (chil (X (n, m) , Y (n, m) )) ; 
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end 


end 

if  atan2 ( Y (n, m) , X (n, m) )  >=  pi/4  &&  atan2  (Y  (n,  m)  ,  X  (n,  m)  )  <  3*pi/4 
if  abs(Y(n,m))  ==  si 

CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI{n,m}  =  zeros(3,3); 


else 

CHI{n,m}  = 

double (chi2 (X(n,m) ,Y (n,m) ) ) ; 

end 

end 

if  atan2 ( Y (n, m) , X (n, m) )  >=  3*pi/4  I  I  atan2  (Y  (n,  m)  ,  X  (n,  m)  )  <  -3*pi/4 
if  abs(X(n,m))  ==  si 

CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI{n,m}  =  zeros(3f3); 


else 

CHI{n,m}  = 

double (chi3 (X(n,m) ,Y (n,m) ) ) ; 

end 

end 

if  atan2 ( Y (n^ m) , X (nf m) )  >=  -3*pi/4  &&  atan2 ( Y (n, m) , X (nf m) )  <  -pi/4 
if  abs(Y(n,m))  ==  si 

CHI{nfm}  =  zeros(3,3); 
elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI{n,m}  =  zeros(3,3); 
elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI{n,m}  =  zeros  (3,3) ; 


else 

CHI{n,m}  = 

double (chi4 (X(n,m) ,Y (n,m) ) ) ; 

end 

end 

CHI_ll(n,m)  =  CHI { nf  m}  (1 , 1 ) ; 
CHI_12(nfm)  =  CHI { n, m}  ( 1 , 2 ) ; 
CHI_22(n,m)  =  CHI { nf m}  (2 f 2 ) ; 
CHI_33(n,m)  =  CHI { n, m}  ( 3 f 3 ) ; 
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G{n,m}  =  besselh(0,l,k_0.*sqrt((X(n,m)-X).A2+(Y(n,m)-Y).A2)); 
G{n,m} (n,m)  =  0; 


end 


end 


i;  0] 

,  [1;  0 

;  0] ) ) ,  [ 

,  [x 

y] ) ; 

[1;  0 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

[1;  0 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

[l;  0 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

i;  0] 

,  [  0 ;  1 

;  0] )  )  ,  [ 

,  [x 

y] ) ; 

[0;  1 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

[0;  1 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

[0;  1 

;  0] ) , 

[X 

y] ) ; 

,  [x 

y] ) ; 

%  Compute  the  derivatives  of  the  susceptibility  appearing  in  (4.60) 
chil_12  =  symfun (con j (dot (chil (x, y) * [0; 1; 0] , [1; i 
chil_12_y  =  symfun (diff (chil_12 (x, y)  ,  y) ,  [x  y] )  ; 
chi2_12  =  symfun (dot (chi2 (x, y) * [0; 1; 0] ,  [1; 0; 0] )  , 
chi2_12_y  =  symfun (diff (chi2_12 (x,  y)  ,  y)  ,  [x  y] )  ; 
chi3_12  =  symfun (dot (chi3 (x, y) * [0; 1; 0] ,  [1; 0; 0] ) , 
chi3_12_y  =  symfun (diff (chi3_12 (x,  y)  ,  y)  ,  [x  y] ) ; 
chi4_12  =  symfun (dot (chi4 (x, y) * [0; 1; 0] ,  [1; 0; 0] ) , 
chi4_12_y  =  symfun (diff (chi4_12 (x, y) , y) , [x  y] ) ; 
chil_22  =  symfun (con j (dot (chil (x,  y) * [ 0;  1 ;  0 ]  ,  [ 0;  : 
chil_22_x  =  symfun (diff (chil_22 (x, y) , x)  ,  [x  y] )  ; 
chi2_22  =  symfun (dot (chi2 (x, y) * [0; 1; 0]  ,  [0; 1; 0] ) , 
chi2_22_x  =  symfun (diff (chi2_22 (x,  y)  ,  x)  ,  [x  y] ) ; 
chi3_22  =  symfun (dot (chi3 (x, y) *[ 0; 1; 0 ],[ 0 ; 1 ; 0 ])  , 
chi3_22_x  =  symfun (diff (chi3_22 (x,  y)  ,  x)  ,  [x  y] )  ; 
chi4_22  =  symfun (dot (chi4 (x,  y)  *[ 0 ;  1 ;  0 ]  , 
chi4_22_x  =  symfun (diff (chi4_22 (x, y)  ,  x) 

CHI_12_y  =  zeros  (size (X,2), size (X, 2)); 

CHI_22_x  =  zeros  (size (X,2), size (X, 2)); 
for  n  =  1 : 1 : size (X, 2 ) 

for  m  =  1 : 1 : size (X, 2 ) 

if  atan2 ( Y (n, m) , X (nf m) )  >=  -pi/4  &&  atan2 ( Y (n, m) , X (nf m) )  <  pi/4 

if  abs(X(n,m))  ==  si 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n/m))  <  si  &&  abs(Y(n/m))  <  si 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n, m)  =  0; 

elseif  abs(X(n/m))  >  s2  |  |  abs(Y(n/m))  >  s2 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n,m)  =  0; 

else 

CHI_12_y  (n,  m)  =  double  ( chil_12_y  (X  (n,  m)  ,  Y  (nf  m)  ))  ; 
CHI_22_x (n, m)  =  double ( chil_22_x (X (n, m) ,  Y (nf m) )) ; 
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165 
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end 


end 

if  atan2 ( Y (n, m) , X (n, m) )  >=  pi/4  &&  atan2  (Y  (n,  m)  ,  X  (n,  m)  )  <  3*pi/4 
if  abs(X(n,m))  ==  si 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI_12_y (n,m)  =  0; 

CHI_22_x (n,m)  =  0; 

else 

CHI_12_y  (n,  m)  =  double  ( chi2_12_y  (X  (n,  m)  ,  Y  (n,  m)  ))  ; 

CHI_22_x (n, m)  =  double ( chi2_22_x (X (n, m) ,  Y (n, m) )) ; 

end 

end 

if  atan2 ( Y (n, m) , X (n, m) )  >=3*pi/4  | |  atan2 ( Y (n, m) , X (n, m) )  <  -3*pi/4 
if  abs(X(n,m))  ==  si 
CHI_12_y (n,m)  =  0; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n,m))  <  si  &&  absfYfn,]!!))  <  si 
CH I_1 2_y ( n , m )  =  0 ; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n,m))  >  s2  |  |  abs(Y(n,m))  >  s2 
CHI_12_y (n, m)  =  0; 

CHI_22_x (n, m)  =  0; 

else 

CHI_12_y (n, m)  =  double ( chi3_12_y (X (n, m) , Y (n, m) )) ; 

CHI_22_x (n, m)  =  double ( chi3_22_x (X (n, m) ,  Y (n, m) )) ; 

end 

end 

if  atan2 ( Y (n, m) , X (n, m) )  >=  -3*pi/4  &&  atan2 ( Y (n, m) , X (nf m) )  <  -pi/4 
if  abs(X(n,m))  ==  si 
CHI_12_y (n,m)  =  0; 

CHI_22_x (n,m)  =  0; 

elseif  abs(X(n,m))  <  si  &&  abs(Y(n,m))  <  si 
CHI_12_y (nfm)  =  0; 

CHI_22_x (nfm)  =  0; 
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elseif  abs(X(n,m) 
CHI_12_y  (n,m) 
CHI_22_x  (n,  m) 

else 

CHI_12_y  (n,m) 
CHI_22_x  (n,  m) 

end 

end 

end 


>  s2  | |  abs (Y (n,m) )  >  s2 
=  0; 

=  0; 

=  double  ( chi4_12_y  (X  (n,  m)  ,  Y  (n,  m) 
=  double  ( chi4_22_x  (X  (n,  m)  ,  Y  (n,  m) 


end 

F  =  cell (size (X,  2 )  ,  size (X,  2 )  )  ; 
for  n  =  1 : 1 : size (X, 2 ) 

for  m  =  1 : 1 : size (X, 2 ) 

F{n,m}  =  pi*li/2*G{n,m}. *exp (li*k_0*X) . * (k_0  A2*CHI_33-li*k_0* . . . 
CHI_22_x+k_0  A2  *CHI_22  +  li*k_0*CHI_12_y )  ; 

end 


end 

El_z  =  zeros (size (X, 2 ) , size (X, 2 ) ) ; 
for  n  =  l:l:size(X,2); 

for  m  =  1 : 1 : size (X, 2 ) 

El_z(n,m)  =  increment A2*trapz (trapz (F { n, m} , 2 )) ; 

end 


end 

E0_z  =  2*pi*exp (li*k_0*X)  ; 

o.  o. 
o  o 


%  Graphs 

surf (X, Y, CHI_1 1 ) 

title (' $\chi_{xx} $  for  Square  Structure interpreter latex ') ; 
xlabel ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 
y label ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 
z label ( ' $\chi_{xx} $ '  ,  ' interpreter ' ,  'latex ' ,  ' font size ',14) ; 
figure 

surf (X, Y, abs (E0_z+El_z) / (2*pi) ) 

title ('1st  Order  Combined  Field  Amplitude  Relative  to  $\left | E_0\right | $ ' , . . . 

' interpreter ' , ' latex ' ) ; 
xlabel ( ' x  (meters ) ' , ' interpreter ' , 'latex ' ) ; 
y label ( ' y  (meters ) ' , ' interpreter ' , 'latex ' ) ; 

zlabel ( ' $\frac{ \left | E_0+E_l\right | } { \left | E_0\right I } $ ' , ' interpreter ' , . . . 
'latex ' , ' font size ' , 14 ) ; 
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221  figure 

222  plot  (X(101,  :  )  ,  abs  (E0_z  (101,  :  )  +El_z  (101,  :  )  )  /  (2*pi)  )  ; 

223  T  =  '1st  Order  Combined  Field  Along  X-axis  Relative  to  $\left | E_0\right | $ ' ; 

224  title  (T,  '  interpreter  '  ,  '  latex  '  )  ; 

225  x  label  (  '  x  (meters )  '  ,  '  interpreter  '  ,  'latex  '  )  ; 

226  ylabel  (  '  $\frac{  \left  |  E_0+E_l\right  |  }  {  \left  |  E_0\right  I  }  $  '  ,  '  interpreter  '  ,  .  .  . 

227  '  latex  '  ,  '  fontsize  '  ,  14  )  ; 


115 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


116 


List  of  References 


[1]  Office  of  Naval  Research.  (2015).  Long  range  broad  agency  announcement 
(BAA)  for  Navy  and  Marine  Corps  science  and  technology.  Internet 
draft.  [Online].  Available:  http://www.onr.navy.mil/~/media/Files/Funding- 
Announcements/BAA/2016/N00014-16-R-BA001.ashx.  Accessed  April  23,  2016. 

[2]  D.  S.  J.  Pendry  and  D.  Smith,  “Controlling  electromagnetic  fields,”  Science,  vol.  312, 
pp. 1780-1782,  June  2006. 

[3]  D.  Schurig  etal.,  “Metamaterial  electromagnetic  cloak  at  microwave  frequencies,” 
Science,  vol.  314,  pp.  977-979,  November  2006. 

[4]  R.  W.  Boyd,  Nonlinear  Optics.  New  York,  NY:  Academic  Press,  2008. 

[5]  Office  of  Naval  Research.  (2015).  Directed  energy  program.  Public  web  page. 
[Online] .  Available:  http://www.onr.navy.mil/Science-Technology/Departments/ 
Code-35/All-Programs/air-warfare-352/Directed-Energy.aspx.  Accessed  April  29, 
2016. 

[6]  G.  Kulacki.  (2014).  An  authoritative  source  on  china’s  military  space  strategy. 
Internet  draft.  [Online].  Available:  http://www.ucsusa.org/sites/default/files/legacy/ 
assets/documents/nwgs/China-s-Military-Space-Strategy.pdf#l.  Accessed  April  29, 
2016. 

[7]  Office  of  Naval  Research.  (2015).  Counter-directed  energy  program.  Public 
web  page.  [Online].  Available:  http://www.onr.navy.mil/en/Science-Technology/ 
Departments/Code-35/All- Programs/aerospace-research-35 1/Counter-  Directed- 
Energy.aspx.  Accessed  April  29,  2016. 

[8]  J.  Pendry  etal.,  “Design  of  electromagnetic  cloaks  and  concentrators  using  form- 
invariant  coodrinate  transformations  of  maxwell’s  equations,”  Photonics  and 
Nanos-tructures  -  Fundamentals  and  Applications,  vol.  6,  pp.  87-95,  Apr.  2008. 

[9]  T.  P.  U.  Leonhardt,  “Transformation  optics  and  the  geometry  of  light,”  Progress  in 
Optics,  vol.  53,  pp.  1-72,  June  2008. 

[10]  O.  Paul  and  M.  Rahm,  “Covariant  description  of  transformation  optics  in  nonlinear 
media,”  Optics  Express,  vol.  20,  pp.  8982-8997,  Apr.  2012. 

[11]  B.  P.  S.  Cummer  et  al.,  “Full-wave  simulations  of  electromagnetic  cloaking  struc¬ 
tures,”  Physical  Review,  vol.  74,  May  2000. 


117 


[12]  A.  K.  W.  Cai,  U.  Chettiar  and  V.  Shalaev,  “Optical  cloaking  with  metamaterials,” 
Nature  Photonics ,  vol.  1,  pp.  224-227,  April  2007. 

[13]  Y.  F.  Y.  Huang  and  T.  Jiang,  “Electromagnetic  cloaking  by  layered  structure  of 
homogenous  isotropic  materials,”  Optics  Express,  vol.  15,  pp.  11  133-11  141, 
September  2007. 

[14]  D.  R.  Smith  et  ah,  “Composite  medium  with  simultaneously  negative  permeability 
and  permittivity,”  Physical  Review  Letters,  vol.  84,  September  2006. 

[15]  E.  J.  Post,  Formal  Structure  of  Electromagnetics:  General  Covariance  and  Electro¬ 
magnetics.  Mineola,  NY:  Dover  Publications. 

[16]  J.  Luscombe,  “Relativity  and  cosmology,”  unpublished. 

[17]  R.  Courant  and  D.  Hilbert,  Methods  of  Mathematical  Physics,  Volume  1.  New  York, 
NY:  Interscience  Publishers,  1954. 

[18]  J.  J.  Sakurai,  Modern  Quantum  Mechanics.  Reading,  MA:  Addison-Wesley  Publish¬ 
ing  Company,  1994. 


118 


Initial  Distribution  List 


1 .  Defense  Technical  Information  Center 
Ft.  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 


119 


